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);
%i=(0:1:(T-1));
%ai=2*pi*i/T-pi;
%M10x(ai(T))=fft(x(ai(T)));
M10x=1/T*fft(x);
L=0;m=length(M10x);
A0=0.2*3.25/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;
for l=1:L-1
if (I(l+1)-I(l))>=sqrt(T)
a(p)=l;
p=p+1;
end
end
a(p)=L;
mm=1;
for q=1:p
max=abs(M10x(I(a(q))));
bi(q)=I(a(q));
for l=mm:a(q)
if abs(M10x(I(l)))>max
max=abs(M10x(I(l)));
Fmax=I(l);
bi(q)=Fmax;
end
end
mm=a(q)+1;
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
%**********************************************************************************************
y0(t)=x(t);
for k=1:m-1
w(k)=y(k);
end
for k=1:m-1
mw(k,:)=fft(y0);
y0=(y0.*exp(-j*w(k)*t)-mw(k,:)).*exp(j*w(k)*t);
end
%for i=1:(t-1)
% ai(T)=2*pi*i/T-pi;
%end
M20x=1/T*fft(x(t).*x(t));
L=0;m=length(M20x);
A0=0.2/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
w(n)=pi*bi(n)/T-pi/2;
end
colperm(w);
fl=colperm(1);
m=1;
for n=1:p
if n~=fl
w(m)=w(n);m=m+1;
end
end
阅读(857) | 评论(1) | 转发(0) |