svd_tls.m

使用自编函数基于奇异值分解总体最小二乘法（svd-tls）实现AR模型谱估计

%使用自编函数svd-tls实现AR模型谱估计

close all;
clear all;
clc;
%tt=zeros(5,1);
N=128;
n=1:128;
avec=0;
%fc=invsp(:,1)/invsp(1,1);%a1,...ap
a1=sqrt(20);
a2=sqrt(2);
%w=randn(128,1);
%x=a1*sin(2*pi*0.2*n)+a2*sin(2*pi*0.213*n)+w';%产生数据序列
%xxcorr=xcorr(x);%自相关序列
%--------------------
SPy11=0;%20次SVD-TLS
VSPy11=0;
for k=1:20
    
    
w=randn(128,1);
x=a1*sin(2*pi*0.2*n)+a2*sin(2*pi*0.213*n)+w';%产生数据序列
xxcorr=xcorr(x);%自相关序列
   
%-----------交换PE和M的值－－－－－－－－－－－－－－
%-----------M>=pe
pe=31;%为超定矩阵
M=47;%colum
for i=1:M
    for j=1:pe+1
        re(i,j)=xxcorr(pe+i+1-j);%生成自相关矩阵（M*(pe+1))
    end
end
disp('自相关矩阵的大小=');
relength=size(re)%
[u,s,v]=svd(re);%svd分解
%-------h=min(m,n)-----
%h=pe=31取最小的
%--------change M to pe
for i=1:pe
    if s(i,i)/s(1,1)>0.05%归一化
        if i>0;
        temp=i;%有效p
        end
    else
        s(i,i)=0;%生成se
       
    end
end
disp('有效秩的个数')
p=temp
%sp
%-------------------------------------------
sp=zeros(p+1);
for j=1;p
    for i=1:pe+1-p
       
        %sp=sp+(s(j,j)^2)*(v([i:i+p],j))*(v([i:i+p],j))';%sp矩阵生成
        sp=sp+(s(j,j)^2)*(v([i:i+p],j))*(v([i:i+p],j))';%sp矩阵生成
        %disp(sp);
    end
end
%disp(sp);
invsp=inv(sp);%求逆
fc=invsp(:,1)/invsp(1,1);%a1,...ap
%-----------------------------joseph add
%avec=fc+avec;
[P,F]=freqz(1,fc,1024,1);
SPy11=SPy11+P;
%F=SPy11+F;
%VSPy11=VSPy11+abs(fft(fc)).^2;
%------------------------------------
freqz(1,fc,1024,1);%显示图形
%plot(fc);
%title('svd-tls');
legend('20次SVD-TLS估计图');
hold on;
end

figure(6)
freqz(1,avec/20,1024,1);%显示图形
legend('20次SVD-TLS估计的均值');
figure(7)
plot(SPy11/20);
legend('20次SVD-TLS估计的均值');