order1.m

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

T=256;
w1=-2.4;w2=-0.85;w3=1.25;
u1=1.8;u2=0;u3=-0.5;
t1=0:0.0001:0.5;
e1=1*exp(-1*t1);
t=2:(T+1);
s1=2.8*e1(t)+0.2*e1(t-1);
e2=0.8*exp(-0.8*t1);
s2=-2*e2(t)-0.5*e2(t-1);
e3=1+randn(1,T+1);
s3=1.8*e3(t)+0.9*e3(t-1);
e4=randn(1,T+1);
v=1.3*e4(t)+e4(t-1);
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);
M10x=1/T*fft(x);               
L=0;m=length(M10x);
A0=0.2*3.25*1/(T^(1/16));
for k=1:m
    if abs(M10x(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(M10x(I(1)));
 bi(1)=I(1);
for l=1:L-1
    if abs(M10x(I(l)))>max
       max=abs(M10x(I(l)));
       Fmax=I(l);
       bi(n)=Fmax;
    end    
    
    if (I(l+1)-I(l))>=sqrt(T)
        p=p+1;max=abs(M10x(I(l+1)));n=n+1;bi(n)=I(l+1);
    end
end
for n=1:p
     w(n)=(2*pi*bi(n)/T)-pi;
end
aaa=min(abs(w));
m=1;
 for n=1:p
    if abs(w(n))~=aaa
       y(m)=w(n);m=m+1;
    end
 end
%**********************************************************************************************