SVD_TLS,Cadzow 谱估计,周期图法的matlab仿真

SVD_TLS,Cadzow 谱估计,周期图法的matlab仿真

• author:hwb
• version:0.01
• 依然是信号处理作业,细节日后再补

现代谱估计的Cadzow谱估计子ARMA模型matlab实现
顺便加了个周期图估计的仿真用作对比

matlab代码

clear;
clc;

%% 参数设置
sub_fre = 3;
fs=512;
th = 0.05;
%% 仿真数据产生
m=0:(fs-1);
f_sub = linspace(0.2,0.4,sub_fre);
A = ones(sub_fre,1);
xn = zeros(fs,1)';
for i=1:sub_fre
    xn = xn + A(i)*sin(2*pi*f_sub(i)*m); %产生含有噪声的序列xn
end
xn = xn + randn(1,fs);
%% 自相关矩阵生成
Rx=xcorr(xn);
pe=100;                   % pe>=p;
qe=150;                   % qe>=q,(qe-pe)>=(q-p);
M=300;                    % M>=pe;
Re=zeros(M,pe+1);
for i=1:M
    for j=1:(pe+1)
       Re(i,j)=Rx(1,fs+qe+i-j+1); % fs 是因为Rx得出来的是偶对称的自相关序列,取后半段即可
    end
end

%% 求AR的阶数 和前p个sigma
[U,S,V]=svd(Re); % 奇异值分解
%%求p值
p=0;
SS=zeros(1,pe+1)';
for i=1:(pe+1)
    if S(i,i)/S(1,1)>th
       p=p+1;
       if p>0
           SS(p)=S(i,i);
       end
    end
end


%% 求Sp及其逆
Sp=zeros(p+1,p+1);
for j=1:p
    for i=1:(pe+1-p)
      Sp=Sp+SS(p)^2*V(i:i+p,j)*V(i:i+p,j)';
    end
end
inv_Sp=inv(Sp);
%% 总体最小二乘估计
x=zeros(1,p+1)';
for i=1:(p+1)
   x(i)=inv_Sp(i,1)/inv_Sp(1,1);
end

%% 求响应
A=x;
% freqz(1,A,fs,1);
% grid on;
% title('SVDTLS')
%% Cadzow 谱估计
Cxk2 = xcov(xn);
Cxk = zeros(1, fs);
for k=1:fs/2
    Cxk(k) = Cxk2(1,fs+k);  % 提取协方差矩阵
end
Cxk = [Cxk(1)/2,Cxk(2:128)]; % 4.4.6,第一个数据取1/2
nk=zeros(1,p+1);
for k=1:(p+1)
    for i=1:p
        a=A(i)*Cxk(abs(k-i)+1); % pk是偶函数
        nk(i)=a+nk(i);
    end
end
%freqz(nk,A,fs,1);

%% 周期图法
window=boxcar(length(xn)); %矩形窗
nfft=fs;
[Pxx,f_d]=periodogram(xn,window,nfft,fs); %直接法

%% 画图
subplot(2,2,1); %原信号
x = 1:fs;
plot(x, xn);
title("原信号")

subplot(2,2,2); %AR过程功率谱
[H1,w1]=freqz(1,A,fs,1); 
Hf1=abs(H1);  %取幅度值实部
plot(w1,20*log(Hf1))  
title("AR过程功率谱(q=0)")

subplot(2,2,3); %Cadzow 谱估计子
[H2,w2]=freqz(nk,A,fs,1); 
Hf=abs(H2);  %取幅度值实部
plot(w2,20*log(Hf)) 
title("Cadzow 谱估计子")

subplot(2,2,4); %Cadzow 谱估计子
plot(f_d/fs,10*log10(Pxx));
title("周期图法")
total_title = ['谱估计方法对比:频率点=' num2str(sub_fre) ' 阈值=' num2str(th) ];
sgtitle(total_title); % matlab2019
%%

仿真结果

暂时没空优化这好纵轴

结语

冲 有兴趣看现代信号处理笔记的朋友可以点个赞.

参考文献

现代信号处理(第三版) - 张贤达

Matlab 实现经典功率谱分析和估计

最小二乘法和SVD_TLS法