LS估计的AR参数Cadzow谱估计信号功率谱密度

clc;
clear;
close all;

n=[1:128];x=sqrt(20)*sin(2*pi*0.2*n)+sqrt(2)*sin(2*pi*0.213*n)+randn(1,128);r11=0;%对信号x(n)进行采样，采样点数为128个
for i=1:127
    for j=i+1:128
   a=x(j)*x(j-i);r11=r11+a;
    end
   r1(i)=r11/128;r11=0;
end
r12=0;
for i=1:128
   a=x(i)^2;r12=r12+a;
end
r1=[r12/128 r1];%计算x的自函数
r2=r1(61:128);
p=4;
r3=fliplr(r1(61-p:61));
RR=toeplitz(r2,r3);%设一个AR阶数p，构造一个M X p+1维的法方程系数矩阵RR
Y=RR'*RR;
 Yinv=inv(Y);
for i=1:p
    pp(i)=Yinv(i+1,1)/Yinv(1,1);
end%用LS算法求出这p个参数
m=mean(x);
c=r1-m*m;
c1=[c(1)/2,c(2:128)];
nk=zeros(1,p+1);
pa=[1 pp];
for i=1:p+1
    for j=1:i
    a=pa(j)*c1(i-j+1);
    nk(i)=a+nk(i);
    end
   end
Nz=0;Az=0;f=0:0.001:0.5;
for j=0:p
   Nz=nk(j+1)*exp((-1)*2i*pi*f*j)+Nz;
   Az=pa(j+1)*exp((-1)*2i*pi*f*j)+Az;
   end
Px=real(Nz./Az)*2;%用Cadzow谱估计子估计出信号的功率谱密度
plot(f,Px);grid on