order3.m

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

T=8192;
w1=-1;w2=-0.5;w3=0.75;
u1=0.6;u2=1.8;u3=-0.5;
t1=0:0.0001:0.9;
e1=0.84*exp(-0.84*t1);
t=2:(T+1);
s1=2.8*e1(t)+0.2*e1(t-1);
e2=0.85*exp(-0.85*t1);
s2=-2*e2(t)-0.5*e2(t-1);
e3=0.845*exp(-0.845*t1);
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);
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
M30x=1/T*fft(y(t).*y(t).*conj(y(t)));
L=0;m=length(M30x);
A0=0.32*4.798*1/(T^(1/16));
for k=1:m
    if abs(M30x(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(M30x(1));
bi(1)=I(1);
for l=1:L-1
    if abs(M30x(I(l)))>max
       max=abs(M30x(I(l)));
       Fmax=I(l);
       bi(n)=Fmax;
   end      
    if (I(l+1)-I(l))>=sqrt(T)
        p=p+1;max=abs(M30x(I(l+1)));n=n+1;bi(n)=I(l+1);
    end
end
  for n=1:p
     ww(n)=(2*pi*bi(n)/T)-pi;
  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