pm2_wv_doa_test.m

非平稳信号处理中用维格纳时频分布进行两个正弦调频信号的波达方向的估计

clear;
clear all;
format short;
c=3*10.^8;
L=8;
N=256*1;
fm1=0.25;
fm2=0.32;
Km1=0.025;
Km2=0.033;
f0=0.32;
lamta=c/f0;
len=lamta/2;
P1=20*pi/180;
P2=40*pi/180;
snr=10;
Amp=sqrt(2*10^(snr/10));
pm1=1/(2*pi);
pm2=1/(2*pi);
t=1:N;
ts=t/512;
sig1=Amp*exp(j*2*pi*(fm1*t-pm1/(2*pi)*cos(2*pi*ts)));
f1=fm1+pm1*sin(2*pi*ts);
[tfrf1,t,f]=tfrsp(sig1(1,:).'); 
figure(1);
[t,f]=meshgrid(t,f);
contour(t,f,tfrf1);
%subplot(3,1,1);
%plot(ts,f1);
sig2=Amp*exp(j*2*pi*(fm2*t-pm2/(2*pi)*cos(2*pi*ts)));
f2=fm2+pm2*sin(2*pi*ts);
%subplot(3,1,2);
%plot(ts,f2);
s=[sig1;sig2];
  for i=1:N;
  l=1:L;  
  x1=exp(j*2*pi*(len*f1(1,i)*(l-1)*sin(P1))/c);
  x2=exp(j*2*pi*(len*f2(1,i)*(l-1)*sin(P2))/c);
  a1=x1.';
  a2=x2.';
  a=[a1,a2];
  S(:,i)=a*s(:,i);
  end
  
noise=randn(L,N)+j*randn(L,N);
z=S+noise;
[tfrf,t,f]=tfrsp(z(1,:).'); 
figure(2);
[t,f]=meshgrid(t,f);
contour(t,f,tfrf);
rou=round(0.3/0.5*N/2);
[Mm1 Nn1]=max(tfrf(1:rou,:));
ff01=mean(Nn1(N/2-5:N/2-5))./N;
lamta01=c/ff01;
[Mm2 Nn2]=max(tfrf((rou+1):N,:));
ff02=(mean(Nn2(N/2-5:N/2-5))+rou)./N;
lamta02=c/ff02;
rou1=round(0.3/0.5*N);
 for i=1:L 
    zz=z(i,:).';
    [tfr,t,f]=tfrspwv([z(1,:).',zz]);
     y1=max(tfr(1:rou1,(N/2-5:N/2-5)));
     yy1(i,:)=y1;
     y2=max(tfr(rou1+1:N,(N/2-5:N/2-5)));
     yy2(i,:)=y2;
 end;
Rz=(yy1*yy1')/N;
J = zhihuan(L);
[e,v]=eig(Rz);
  if v(1,1)<v(L,L)
     vv=J*v*J;
     ee=J*e*J;
  else  vv=v;
     ee=e;
  end
es=ee(:,1:1);
en=ee(:,(2:L));
numm=900;
aaaa=zeros(L,numm);
 for k=1:L;
     for h=1:numm;
         aaaa(k,h)=exp(j*2*pi*(k-1)*len*ff01*sin((0.1*h*pi/180))/c);
     end
 end
b=eye(L);
cc=en*en';
 for m=1:numm;
     aac=aaaa(:,m);
     pp=aac'*aac/((aac)'*cc*aac); 
     p1(m)=real(pp);
 end
 [m1,n1]=max(p1);
 doa1=n1/10
Rz=(yy2*yy2')/N;
[e,v]=eig(Rz);
J = zhihuan(L);
  if v(1,1)<v(L,L)
     vv=J*v*J;
     ee=J*e*J;
  else  vv=v;
     ee=e;
 end
 es=ee(:,1:1);
 en=ee(:,(2:L));
aaaa=zeros(L,numm);
 for k=1:L;
     for h=1:numm;
         aaaa(k,h)=exp(j*2*pi*(k-1)*len*ff02*sin((0.1*h*pi/180))/c);
     end
 end
b=eye(L);
 %cc=b-es*es';
 cc=en*en';
 for m=1:numm;
     aac=aaaa(:,m);
     pp=aac'*aac/((aac)'*cc*aac); 
     p2(m)=real(pp);
 end

[m2,n2]=max(p2);
doa2=n2/10
u=0:0.1:90-0.1;
%subplot(3,1,3);
figure(3);
plot(u,p1,u,p2);
xlabel('角度（度）');
ylabel('功率谱');