谐波恢复的ARMA建模算法的Matlab实现

谐波恢复的ARMA建模算法的Matlab实现

源自:Pisarenko 谐波分解法

• 依然是信号处理作业,细节日后再补

谐波恢复的ARMA建模算法

% 谐波恢复的ARMA建模算法
clear;
clc;
%% 参数设置
sub_fre = 2;
fs=512;
th = 0.1;
%% 仿真数据产生
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=10;                   % pe>=2p;
M=40;                    % M>>p;
Re=zeros(M,pe+1);
for i=1:M
    for j=1:(pe+1)
       Re(i,j)=Rx(1,fs+pe-j+i+1); % Rx是偶对称的自相关序列,取后半段
    end
end
%% 利用SVD_TLS算法得到2p 和a
[A,p] = func_SVD_TLS(Re, pe, th);
a_sym = A;
p = p/2;
%% 计算特征多项式
z_res = roots(a_sym);
%% 计算各个谐波频率
f_res = []; 
for i=1:2*p % 有些时候不太稳定,只能全跑一遍了.
    f_course = atan(imag(z_res(i))/real(z_res(i)))/pi;
    if(f_course > 0)
        f_res = [f_res, f_course];
    end
end
f_res = sort(f_res);
disp(f_res);
%%

SVD_TLS整体最小二乘估计

其中,SVD_TLS部分被封装成函数:

function [A,p] = func_SVD_TLS(Re, pe, th)
%FUNC_SVD_TLS 此处显示有关此函数的摘要
% pe>=p; % qe>=q,(qe-pe)>=(q-p);% M>=pe;
%% 求AR的阶数 和前p个sigma
[~,S,V]=svd(Re); % 奇异值分解
%%求p值
p=0;
SS=zeros(1,pe+1)';
for i=1:pe
    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=pinv(Sp);
%% 总体最小二乘估计
x=zeros(1,p+1)';
for i=1:(p+1)
   x(i)=inv_Sp(i,1)/inv_Sp(1,1);
end
%% 求响应
A=x;
end

详情可见:

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

结语

匆忙得没空写结语

参考文献

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