chrip1_music_rmusic.m

基于MUSIC算法的doa估计算法的改进,比原来的算法更有效,大家可以下载看一下了,欢迎讨论.

  %  the  data of signal
clear;
format long;
c=3*10.^8;
L=9;
sam=128*4;
N=128*4;
f0=0.175;
lamta=c/f0;
len=lamta/2;

P1=20*pi/180;
doature=20;
for  nn=1:6;   
  snr=5*(nn-2);
  Amp=sqrt(2*10^(snr/10));
%snr=10;
%Amp=sqrt(2*10^(snr/10));
sig=Amp.*fmlin(N,0.1,0.25);
s=sig.';
tt=1:N;
f1=(0.15/N)*tt+0.1;

S=zeros(L,N);
 
  for t=1:N;
  i=1:L;  
  x1=exp(j*2*pi*(len*f1(1,t)*(i-1)*sin(P1))/c);
  a1=x1.';
  a=a1;
  S(:,t)=a*s(:,t);
  end
  for nnn=1:10; 
      noise=randn(L,N)+j*randn(L,N);
      z=S+noise;
 
 %计算中心频率
      [tfrf,t,f]=tfrsp(z(1,:).');
      [Mm Nn]=max(tfrf);
      ff0=mean(Nn(N/2-5:N/2-5))./N;
      lamta0=c/ff0;
      %f01=min(Nn)./N;
      %f02=max(Nn)./N;
      %f0=(f01+f02)./2;
 %计算阵元间的互WV分布
      for i=1:L 
          zz=z(i,:).';
          [tfr,t,f]=tfrwv([z(1,:).',zz]);
          %[WH,rho,theta]=htl(tfr,N,N);
          %figure(i);
          %contour(t,f,tfr);
          %y=max(tfr(:,N/2-5:N/2+5));
          y=max(tfr);
          %y=max(WH);
          yy(i,:)=y;
      end;
          %mesh(t,f,abs(tfr));
          Rz=(yy*yy')/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,300);

          for k=1:L;
              for h=1:300;
                  aaaa(k,h)=exp(j*2*pi*(k-1)*len*ff0*sin((0.1*h*pi/180))/c);
                 
             end
          end
          b=eye(L);
          cc=b-es*es';
          for m=1:300;
              aac=aaaa(:,m);
              pp=aac'*aac/((aac)'*cc*aac); 
              p(m)=real(pp);
          end
 
 %p1=p(1:450);
 %p2=p(451:900);
%[m1,n1]=max(p1);
%doa1=n1/10
         [m,n]=max(p);
         doa(nnn)=n/10;
         
         
           %构造多项式
          Un1=ee(1,2:L);
          Un2=ee(2:L,2:L); 
          Li=[1 0 0 0 0 0 0 0];
         c1=Un1*inv(Un2)*Li.';
         ccc=[1 c1];
         R=roots(ccc);
         est1=asin(angle(R)*lamta0/(2*pi*len));
         agl1(nnn)=-1*est1*180/pi;
  
      end
aav=mean(doa)
bbv=mean(agl1)
%doaa=doa-doature;
%aasss(nn)=sqrt(sum(doaa.*doaa)./nnn);
aas=std(doa);
bbs=std(agl1);
aavv(nn)=aav;
aass(nn)=aas;
bbss(nn)=bbs;

end
hold on;
i=1:6;
snr=(i-2)*5;
%semilogy(snr,aass,'b-*',snr,aasss,'r-*');
semilogy(snr,aass,'b-*',snr,bbss,'r-*');
%text(21,0.2,'70度');
%text(21,0.05,'30度');
xlabel('信噪比SNR(dB)');
%ylabel('估计角度的均方根误差');
legend('标准差','均方根误差RMSE')