tops.m

宽带信号源测向的一种新方法,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;