1_2_3_xiebo.rar_harmonic order

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