fangzhen1_ls.m

此程序为MATLAB实现的现代信号处理中几种常用的谐波恢复方法

function fangzhen1_ls()

pe=6;
%生成增广矩阵Re
M=100;
loops=20;
%pe=25;
%Re=zeros(M,pe+1);
fvMatrix=zeros(pe/2,loops);
arMatrix=zeros(pe,loops);
for loop=1:1:loops
%    w = wgn(1,2000,0);
    w = randn(2000,1);
    for n=1:1:128
        x(n)=(20^(1/2))*sin(2*pi*0.2*n)+(2^(1/2))*sin(2*pi*0.213*n)+w(n+50*loop);
        %x(n)=(20^(1/2))*sin(2*pi*0.2*n);
    end
    Rxx=xcorr(x,'unbiased');
    for i=1:1:M
        for j=1:1:pe+1
            Re(i,j)=Rxx(pe+i+1-j+128);
        end
    end

    kk=pe;
    %最小二乘法估计ARMA模型的参数
    b=-1*Re(:,[1]);
    A=Re(:,[2:1:kk+1]);
    %ar_lsline=(A'*A)\(A'*b);
    ar_lsline=inv(A'*A)*A'*b;
    arMatrix(:,loop)=ar_lsline;
    
    %求谐波频率
    for j=2:1:kk+1
        fc(j) = ar_lsline(j-1);
    end
    fc(1) = 1;
    fz=roots(fc);
    count=1;
    for j=1:2:kk
        fw(count)=atan(imag(fz(j))/real((fz(j))));
        fv(count)=fw(count)/(2*pi);
        count=count+1;
    end
    fv=abs(fv);
    fv=sort(fv);
    fvMatrix(:,loop)=fv';
end

%分别求AR参数估计和频率估计的均值和方差
arMatrix
ar_mean = mean(arMatrix')'
ar_var = var(arMatrix')'
fvMatrix
fv_mean = mean(fvMatrix')'
fv_var = var(fvMatrix')'