Chinaunix首页 | 论坛 | 博客
  • 博客访问: 303268
  • 博文数量: 71
  • 博客积分: 2010
  • 博客等级: 大尉
  • 技术积分: 734
  • 用 户 组: 普通用户
  • 注册时间: 2009-05-20 17:17
文章分类

全部博文(71)

文章存档

2011年(1)

2009年(70)

我的朋友

分类:

2009-06-10 16:30:55

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) |
0

上一篇:MatlabTest060902

下一篇:DTFT

给主人留下些什么吧!~~

chinaunix网友2009-09-07 23:25:22

没注释啊,算什么的