参考论文：《相干信源波达方向估计的加权空间平滑算法》

王布宏 王永良 陈辉

武汉空军雷达学院重点实验室 武汉430010

摘  要 提出了一种用于空间相干源 DOA 估计的加权空间平滑算法 WSS, weighted spatial smoothing

常规的空间平滑算法没有利用子阵输出的互相关信息 而且对相干信源的分辨力较差 WSS 算法充分

利用了子阵输出的自相关信息和互相关信息 将主阵协方差矩阵的所有子阵阵元数阶的子矩阵进行加

权平均 以期提高常规空间平滑算法的分辨性能 文中以加权平滑后等效的信源协方差矩阵的对角化

为约束条件 推导了加权矩阵的理论表达式 计算机仿真结果表明 WSS 算法与常规空间平滑算法相

比具有更高的分辨性能和更低的信噪比门限 特别是在子阵划分较多时其优越性更加明显

关于这篇论文的算法我做了一个仿真，输入两个相干信源，如果两个信源角度相邻(<30，还有许多很糟的情况),算法仿真完全失效。

有人在这篇论文有过研究没？不知道是我的代码问题，还是这篇论文数据有问题(本人水平低，这种情况应该不可能)，贴出代码与论文，大家帮我分析一下：

%DOA估计---加权空间平滑music算法,使用加权前后向空间平滑clear all;close all;clc;

derad = pi/180;        % deg -> radradeg = 180/pi;        %rad—>degsource_number=3;  %信源数sensor_number=8 ; %阵元数N_x=200;        %信号长度snapshot_number=N_x;%快拍数/阵元采样数f0=1e6;w0=f0*2*pi;fs=10*f0;ts=1/fs;ratio_dl=0.5;  %阵元间距与波长的比值l=0.15;  %中心频率对应的信号波长，无所谓d=ratio_dl*l;  %阵元间距：必须小于1/2波长snr=20;   %信噪比dd=0:d:(sensor_number-1)*d;        % 阵元相对位置矢量

source_theta=[10 35 40];   %信号的入射角% s0=randn(1,N_x);% s=[s0;s0];s=[exp((w0*1i*[0:N_x-1])*ts);exp((w0*1i*[0:N_x-1])*ts);exp((w0*1i*[0:N_x-1])*ts)];      %仿真相干的入射信号A=exp(-1i*2*pi*dd.'*sin(source_theta*derad)/l);              %生成导向矩阵x=A*s;     %阵元接收信号X=awgn(x,snr,'measured');   %阵元接收信号时加高斯白噪声

Rx=X*X'/snapshot_number;   %波束形成器接收信号的协方差矩阵

%----------------------------(FBSS)前后向空间平滑算法、(直接由协方差矩阵进行平滑)----------------------------------p1=5;Rxx_s= fb_spatical_smoothing(Rx, p1);   %空间平滑处理对Rxx降维[U0,S0,V0]=svd(Rxx_s);  %特征值分解 Rxx = U*S*V' 其中特征值按大小降序排列在对角线上Un0=U0(:,source_number+1:p1); %取出小特征值对应的特征向量构成噪声子空间Gn0=Un0*Un0';

searching_theta=-90:0.1:90;%线阵的搜索范围为-90~90度 for i=1:length(searching_theta)   a_theta0=exp(-1i*2*pi*dd(1:p1)'*sin(searching_theta(i)*derad)/l); %搜索矢量，列向量8*1   Power_temp0(i)=abs((a_theta0)'*Gn0*a_theta0);   Pmusic0(i)=1/Power_temp0(i); end normPower_out0=10*log(Pmusic0/max(Pmusic0));  %对输出频谱归一化处理%---------------------------------------------------------------------------------------------------

%------------------------------(WFBSS)加权前后向空间平滑----------------------------------------------------------%--1.对协方差矩阵进行信号预加强处理-----------------------------------------------------------------[V0,D0]=eig(Rx);    %特征值升序排列D00=diag(D0);   %取出特征值Dn=D00(1:sensor_number-1);%取出噪声对应的特征值这里相干信源组数为1Dnm=mean(Dn);Dm=D00(sensor_number)-Dnm;Rxm=V0(:,sensor_number)*Dm*V0(:,sensor_number)';

%--2.CBF预估计角度-----------------------------------------------------------------------------------searching_doa=-90:0.1:90;  %对线阵的角度搜索范围为-90~90度，精度为0.01度 for i=1:length(searching_doa)   Wa=[exp(1i*dd*2*pi*sin(searching_doa(i)*derad)/l)].';%每个扫描角度对应的权向量   CBF_Power_out(i)=abs(Wa'*Rxm*Wa); end [pks,locs]=findpeaks(CBF_Power_out);   %找极大值与极大值坐标，不是谱峰位置，包括功率谱中其他峰值点mean_pks=mean(pks);  %求出峰值的平均值m_pks=find(pks>mean_pks);%找出pks中较大的值得位置n_m_pks=length(m_pks);%pks中较大值的个数，也就是doa_locs=zeros(1,n_m_pks);for i_pks=1:n_m_pks    doa_locs(i_pks)=locs(m_pks(i_pks));end% CBF_doa=(doa_locs-1)*0.1-90;%CBF粗估计出的角度 CBF_doa=source_theta;   %仿真直接取真实角度值%--3.构造矩阵B----------------------------------------------------------------------------------------p2=5; %子阵阵元数K=sensor_number+1-p2;%子阵数B=exp(-1i*2*pi*dd(1:K).'*sin(CBF_doa*derad)/l);%K*M，K子阵数，M信源数%--4.利用伪逆求最优权矩阵a0=0.2;aa=a0.^([0:K-1]);    %向量% W=toeplitz(aa);% W=pinv(B*B');W=eye(K);%--5.加权空间平滑Rwfb=W_fbss(Rx,W);%--6.music算法求谱峰[V1,D1]=eig(Rwfb);  %特征值分解 其中特征值升序排列在对角线上Un=V1(:,1:p2-source_number); %取出小特征值对应的特征向量构成噪声子空间Gn=Un*Un';

searching_theta=-90:0.1:90;%线阵的搜索范围为-90~90度 for i=1:length(searching_theta)   a_theta=exp(-1i*2*pi*dd(1:p2)'*sin(searching_theta(i)*derad)/l); %搜索矢量，列向量4*1   Power_temp(i)=abs((a_theta)'*Gn*a_theta);   Pmusic(i)=1/Power_temp(i); endnorm_Power_out=10*log(Pmusic/max(Pmusic));  %对输出频谱归一化处理

plot(searching_doa,normPower_out0,'b',searching_doa,norm_Power_out,'r');%axis([-90 90 -90 90]);xlabel('入射角/degree');ylabel('归一化空间谱/dB');legend('前后向空间平滑','加权前后向空间平滑');title('加权空间平滑DOA估计(多信源输入)');grid on;

2014-7-20 16:34 上传

点击文件名下载附件

116.95 KB, 下载次数: 8656