TOPS，宽带信号源测向的一种新方法

% --------------------------------------------------------------------------------------------------- % TOPS method % --------------------------------------------------------------------------------------------------- clear all Nsensor=12; Radius=0.5; %注意这里Radius仅仅是与中心频率f0对应波长的比值 DOA=[10 14 40]; Nsignal=length(DOA); SNR=13*[1 1 1]; snapshot=20; fl=200; f0=300; fh=400; Fs=1000; bw=fh-fl; Nfft=128;%选取每段的时域数据长度，即dft长度 N=Nfft*snapshot;%产生时域信号的总点数 Pn=0.01;%因为噪声肯定一样，所以要令其固定，来确定不同信号源的功率 Ps=Pn*10 .^(SNR/10); % f_focus=[]; % X_focus=[]; % result=0; %下面产生高斯白噪声，并且低通滤波得到宽带信号，存放在st中 ss=randn(N,Nsignal)*diag(sqrt(Ps)); st=[]; [bb,aa]=butter(10,[fl/(Fs/2),fh/(Fs/2)]); for r=1:Nsignal temp=filter(bb,aa,ss(:,r)); st=[st,temp]; end %产生信号的功率谱 window=[];%求功率谱加的窗 noverlap=[]; [Pxx,f]=psd(ss(:,1),Nfft,Fs,window,noverlap); figure,plot(f,Pxx); [Pyy,f]=psd(st(:,1),Nfft,Fs,window,noverlap); figure,plot(f,Pyy); J1=round(Nfft*fl/Fs+1);%最低频率fl对应DFT的第J1个点 J2=round(Nfft*fh/Fs+1); J=J2-J1+1;%选用的频率点数目 fll=(J1-1)*Fs/Nfft;%第J1个点对应的模拟频率 fhh=(J2-1)*Fs/Nfft;%第J2个点对应的模拟频率 subband=linspace(fll,fhh,J);%各个子带频点 Xfallsnapshot=zeros(Nsensor,J,snapshot);%用于存放多次快拍频域数据，Xfallsnapshot为3维矩阵 for number=1:snapshot%快拍循环开始 Sf=fft(st((number-1)*Nfft+1:number*Nfft,:)); %fft是按列进行的，每一列都变换，对应各个信号 %Sf的规模为（Nfft，Nsignal） Sf=Sf(J1:J2,:);%将Sf频率分量为零的部分去掉，保留规模为（J，Nsignal） %产生出复噪声 randn('state',sum(100*clock)); nn=sqrt(Pn/2)*randn(Nsensor,J); randn('state',sum(100*clock)); nn=nn+j*sqrt(Pn/2)*randn(Nsensor,J); % nn=randn(Nsensor,J); %下面是把白噪声变成色噪声，仍然存在nn中 % nt=[]; % for rr=1:Nsensor % temp=filter(bb,aa,nn(rr,:));%滤波函数filter对于行或者列数据都可以，且还保持原来的规模，即行仍为行，列仍为列 % nt=[nt;temp]; % end % nn=nt; Nf=fft(conj(nn'));%因为fft对列进行，而这里需要对nt的各行进行fft，所以先对nt进行转置，fft后再转置回来 Nf=conj(Nf');%Nf的规模为(Nsensor,J) %下面建立模型得到频域接受数据表达式，存放在Afsf中 Afsf=zeros(Nsensor,J); for m=1:Nsignal a=steerMultiFre(Nsensor,Radius,DOA(m),f0,subband);%a为规模为（Nsensor，J） Sf1=repmat(conj(Sf(:,m))',Nsensor,1);%将每一个信号变成(Nsensor,J)的规模 as1=a.*Sf1; Afsf=Afsf+as1; %各个信号累加得到频域信号数据 end % Xf=Afsf; Xf=Afsf+Nf;%接收到的频域数据 Xfallsnapshot(:,:,number)=Xf; end%快拍循环完毕 %下面是用TOPS法进行DOA估计 X0=Xfallsnapshot(:,(J/2),:); X0=squeeze(X0);%此函数将维数为1的去掉，得到单一频率点(近似f0)处多次快拍频域接收数据矩阵 [W0,F0]=Dmatrix(X0,Nsensor,Nsignal);%对多次快拍用music方法，形参X0为矩阵，每一列为一次快拍，返回噪声和信号子空间 theta=-90:0.1:90; result=zeros(1,length(theta)); for m=1:length(theta) phi=theta(m); D=[]; for n=1:J-1 X=Xfallsnapshot(:,n,:); X=squeeze(X);%此函数将维数为1的去掉，得到单一频率点的多次快拍频域接收数据矩阵 ff=subband(n);%进行估计的频率点 [W,F]=Dmatrix(X,Nsensor,Nsignal); Psensor=(0:1:Nsensor-1)*Radius; Q=exp(-j*2*pi*Psensor*(ff-f0)/f0*sin(phi*pi/180)); Q=diag(Q); U=Q*F0; a3=exp(-j*2*pi*Psensor'*ff/f0*sin(phi*pi/180)); P=eye(Nsensor)-1/(a3'*a3)*a3*a3'; % P=eye(Nsensor); U=P*U; UW=U'*W; D=[D,UW]; end [UU,SS,VV] = svd(D); SS=diag(SS); result(m)=1/min(SS); %取最小的奇异值的倒数来度量是不是真实DOA方向 end result=result/max(result);%归一化 figure, plot(theta,result,'k'); grid on; plotxline(DOA,'m','--') % result=10*log10(result); % figure, % plot(theta,result,'k'); % grid on;