gonglvpu1.m

输入信号为谐波信号

pp=zeros(50,101);
for n=1:50
    
  
A1=1;
A2=4;
ap(1)=1;
N=1001;
fs=N-1;
f1=100;
f2=400;
t=0:1/fs:1;
x=A1*sin(2*pi*f1*t)+A2*sin(2*pi*f2*t)+randn(size(t));
e1=x;
e2=x;
D=2*sumsqr(x)-x(1).^2-x(N).^2;
p=32;
for J=2:p
      A=e1(:,J:N)*e2(:,J-1:N-1)';
      gama=-2/D*A;
      sse=D*(1-abs(gama).^2);
      
    if sse<=0
        break,
    end
        a=e1;
        b=e2;
        e1(:,J:N)=a(:,J:N)+gama*b(:,J-1:N-1);
        e2(:,J:N)=b(:,J-1:N-1)+conj(gama)*a(:,J:N);
        D=D*(1-abs(gama).^2)-abs(e1(J)).^2-abs(e2(N)).^2;
        ap=[ap,0]+gama*conj([0,fliplr(ap)]);
    
 end
 k=0:100;
 W=pi/100*k;
 fenmu=ones(1,101);
 for i=1:J
     fenmu=fenmu+ap(i)*exp(-j*W*(i-1));
 end
 fenzi=sqrt(sse);
 powerf=(abs(fenzi./fenmu)).^2;
 grid on;
 hold on
plot(W/pi,10*log10(powerf));
clear
clc
end