order2.m

原创：在原有的基础上又变成了1_3阶谐波估计

T=1024;
w1=-1;w2=0.5;w3=1;
u1=0.6;u2=1.8;u3=-0.5;
%t1=0:0.0001:0.5;
e1=1+randn(1,T+1);
%e1=1*exp(-1*t1);
t=2:(T+1);
s1=1*e1(t)+0.8*e1(t-1);
e2=1+randn(1,T+1);
%e2=0.8*exp(-0.8*t1);
s2=0.8*e2(t)+0.6*e2(t-1);
e3=1+randn(1,T+1);
%e3=1*exp(-1*t1);
s3=1*e3(t)+0.5*e3(t-1);
e4=randn(1,T+2);
t=3:(T+2);
v=0.96*e4(t)-0.8*e4(t-1)+0.16*e4(t-2);
t=1:T;
x=s1(t).*exp((w1*t+u1)*j)+s2(t).*exp((w2*t+u2)*j)+s3(t).*exp((w3*t+u3)*j)+v(t);
m=3;
w(1)=w1;w(2)=w2;w(3)=w3;
y(t)=x(t);
for k=1:m
    mwk(k)=(1/T)*sum(y(t).*exp(-j*w(k)*t));    
    y(t)=(y(t).*exp(-j*w(k)*t)-mwk(k)).*exp(j*w(k)*t);
end
M20x=1/T*fft(y(t).*y(t));
L=0;m=length(M20x);
A0=0.25*1/T*sum(y(t).*conj(y(t)))*1/(T^(1/16));
%A0=0.3*3.5625*1/(T^(1/16));
for k=1:m
    if abs(M20x(k))>A0
        L=L+1;
        if k<=T/2
            I(L)=k+T/2;
        else
            I(L)=k-T/2;
        end
    end
end
I=sort(I);
p=1;n=1;
max=abs(M20x(1));
bi(1)=I(1);
for l=1:L-1
    if abs(M20x(I(l)))>max
       max=abs(M20x(I(l)));
       Fmax=I(l);
       bi(n)=Fmax;
   end      
    if (I(l+1)-I(l))>=sqrt(T)
        p=p+1;max=abs(M20x(I(l+1)));n=n+1;bi(n)=I(l+1);
    end
end
  for n=1:p
     ww(n)=pi*bi(n)/T-pi/2;
  end
aaa=min(abs(ww));
m=1;
 for n=1:p
    if abs(ww(n))~=aaa
       www(m)=ww(n);m=m+1;
    end
 end