tct2.m

此源码实现DOA估计中的TCT算法

clear all; 
clc;



t=1/400:1/400:40.96; 
f0=100;
fb=40;
tt=0.01;
s1=10*(sin(2*pi*fb*t)./(pi*t)).*exp(j*2*pi*f0*t);                   %宽带信号1为中心频率100hz，带宽40hz的信号    
s2=10*(sin(2*pi*fb*(t+tt))./(pi*(t+tt))).*exp(j*2*pi*f0*(t+tt));    %宽带信号2为宽带信号1延时1/500s
                              
angle1=20*pi/180;
 Ale=20;           
angle2=45*pi/180;
M=1:8;                                  
 
c=300000000;                            
l=c/f0;                                   
d=l/2;                             

for i=1:64
     ss1(i,:)=s1((1+256*(i-1)):(256+256*(i-1)));          
     ss2(i,:)=s2((1+256*(i-1)):(256+256*(i-1)));           %总测时间为40.96s，分成64个段，每段观察时间为0.64s，每段有256个点
end

f00=100;                                   %选取参考频率点

 aa1=exp(j*(M-1)*2*pi*f00*d*cos(angle1)/c);
 aa2=exp(j*(M-1)*2*pi*f00*d*cos(angle2)/c);
 Aa=[aa1.',aa2.'];                                         %参考频率所对应的信号的方向矢量
SNR=15;

  
      for ij=1:64
          a1(ij,:)=exp(j*(M-1)*2*pi*ij*d*cos(angle1)/c);
          a2(ij,:)=exp(j*(M-1)*2*pi*ij*d*cos(angle2)/c);
          A=[a1(ij,:).',a2(ij,:).'];
          N=0.0001*sqrt(1/(10^(SNR/10)))*randn(length(M),length(s1)/64);        
          s=[ss1(ij,:).',ss2(ij,:).'];
          S=s.';
          x=A*S+N;                                             %接收到的信号
      end
for m=1:8
        X(m,:)=fft(x(m,:),64);                                 %将信号进行64点的fft变换
   end
  for ij=1:64                                                  %频带分成64个频率分量
    a1(ij,:)=exp(j*(M-1)*2*pi*ij*d*cos(angle1)/c);
    a2(ij,:)=exp(j*(M-1)*2*pi*ij*d*cos(angle2)/c);
    AA=[a1(ij,:).',a2(ij,:).'];
  end
R11=AA*S*S'*AA';
R12=Aa*S*S'*Aa'
 
 [U1,ss,V1]=svd(R11);
 [U2,ss,V2]=svd(R12);
 T=U2*U1';              %获得聚焦矩阵T
 Y=T*X;            %对接收信号进行聚焦，聚焦到中心频率上
  

Yy=zeros(8,64);
for jji=1:8
        Yy(jji,:)=Yy(jji,:)+Y(jji,:);   
end
YYy(:,:)=Yy(:,:)/64;
      
 YY=YYy(:,1:10);
 Rx=YY*YY'/length(s1);

[V,D]=eig(Rx);
for i=1:length(M)
    u(i)=D(i,i);
end

[B,IX]=sort(u);
for i=1:length(M)-2
    En(:,i)=V(:,IX(:,i));
end

am=exp(j*(M-1)'*2*pi*d*cos(pi*[0.1:0.1:180]/180)/l);
p=En'*am;
for i=1:1800
    p1=norm(p(:,i));
    P(i)=1/p1^2;
end
t=linspace(0.1,180,1800);
PP=10*log10(P);
Maa=zeros(1,1800)          
            
            for ii=1:1798
                if (PP(ii)<PP(ii+1)) && (PP(ii+1)>PP(ii+2))&&(PP(ii+1)>0)
                      
                       Maa(1,ii+1)=PP(ii+1);
                       indices = find(Maa)/10;
                    
                end  
            end
  figure(2);         

 plot(t,10*log10(P));
 
xlabel('波达角/度');
ylabel('空间谱/dB');
  
 e=indices(1)-Ale;
 E(SNR)=sqrt(e*e');


EE=E(SNR1);
 figure(3);
 plot(SNR1,EE); axis([1,15,-1,3]);
xlabel('SNR(dB)');
ylabel('RMSE');
title('RMSE in degrees versus SNR');