随机信号的 DOA 估计方法

1、 随机信号的DOA估计方法一、 实验目的1、 掌握利用周期图法、Capon 方法、MUSIC方法实现随机信号DOA估计的方法。2、 观察阵元数目、阵元间距、信噪比、入射方向等参数对角度谱估计性能的影响。3、 理解特征结构类方法进行DOA估计的优点。二、 实验原理信号的来波方向（DOA）估计石阵列信号处理领域中的一个重要内容。阵列信号模型如图1所示，设均匀阵列中有M个阵元，阵元间距为D，记信号波长为，则阵列等效孔径为，表示入射到阵列上信号的来波方向（DOA），以信号传播方向与阵列法线方向的夹角来表示（顺时针方向为正）。图 1阵列信号模型假设信号源位于远场，即信号在到达各个阵元时的波前为平面波，以

2、原点处的阵元为参考点，则个阵元接收到的信号为 （1）其中，为信号中心频率，为波长。对于窄带解析复信号，有（2）其中为角频率，则第i个阵元上收到的信号可以表示为（3）如果有d个入射源信号，它们的入射角分别为，则有（4）M个阵元接收到的信号用矩阵表示为（5）其中将矩阵写成如下形式，这里为导向矢量。信号的DOA估计大多采用搜索夫人方法，通过对谱估计函数进行峰值搜索，估计信号波数到达的方向。本实验将对周期图法、Capon法以及MUSIC方法予以讨论。1 周期图法已知接收信号观测样本序列为有限长序列，记，其自相关矩阵为。本次实验中根据各态历经假设，对次快拍求平均估计自相关矩阵，从而有。使用周期图方法进行

3、角度谱估计的结果为（6）（7）因此可以通过谱峰搜索估计信号的波达方向。2 Capon方法Capon方法是一种利用空域滤波实现DOA估计的方法，通过在期望方向形成波束，并利用剩余自由度在干扰方向形成零陷从而一致干扰和噪声。为了保证输出信号中包含期望方向的信号，通常保持其外那个方向的增益为1，在此基础上使输出信号的功率最小，输出功率反映了信号的空间谱特性。Capon方法可以描述成一个优化问题：在约束条件或的前提下，使输出功率最小，采用拉各朗日乘子法可得最优解为，对应的角度谱估计结果为（8）其中。3 MUSIC方法（5）式的接收信号形式中，为入射信号，为白噪声，如果个新号批次独立，且与噪声不相关，则

4、有（9）其中。自相关矩阵的特征值和特征向量为：（10）因此可将表示为：，将特征值按从大到小排序，即：其特征结构性质为：（1）（11）（2）（12）个大特征值对应的特征向量被称为信号特征向量，而个小特征值对应的特征向量被称为噪声特征向量，它们张成的空间是信号特征向量张成的空间的正交补。MUSIC方法将接收信号的自相关矩阵进行特征分解，得到信号子空间和噪声子空间，利用噪声特征向量和信号特征向量的正交关系，通过谱峰搜索可以估计到达波方向，从而得到谱估计结果为（13）其中是由噪声特征向量张成的噪声子空间。搜索寻找谱峰来估计到达角度。三、 实验步骤1、仔细阅读有关角度谱估计的内容，按照图2给出的流程图编

5、写使用周期图法进行角度谱估计程序。选择参数M=8,d=2=0.5,N=100,SNR=10dB,中心频率f=300MHz,两个信号源的入射角度分别为-20度和30度，观察谱估计图如图1所示。图1：周期图角度谱估计如图1所示，谱峰位于-20°和30°处，与信号源的入射角度相符。2、观察并记录参数变化对角度谱估计性能的影响。（1） 改变阵元间距d=4=0.25, d=2=2，其他输入同步骤 1，观察谱峰位置是否正确，以及阵元间距对谱分辨率的影响。d=4 d=2图2：周期图角度谱估计（阵元间距）对于空间采样，与奈奎斯特采样定理类似，必须满足d2，所以d=4时，谱峰位置正确，而d=

6、2时，类似欠采样，除了在正确的角度上出现谱峰外，在其他多处地方也出现了谱峰，所以谱峰位置不正确。而由d=4与d=2时的谱峰对比可知，d=2时分辨率更高，因此在满足d2的条件下，d越大，谱分辨率越高。（2）当只有一个入射源信号时，改变=0°,45°,90°，其他输入同步骤1，观察入射方向对角谱估计性能的影响。=0° =45°=90°图3：周期图角度谱估计（入射角度）由图3可知，在=0°,45°时，角度谱上在相应位置出现正确的谱峰，而在=90°时，在90°和-90°处都出现一个谱峰，影响D

7、OA估计心能。不过，在实际中，90°和-90°物理朝向相差明显，可以人工判断到底是90°还是-90°。（3）改变M=4，16，其它条件同步骤1，观察阵元数目对谱估计性能的影响。M=4 M=16图4：周期图角度谱估计（阵元数目）由图4可知，M=4时，相对于M=8的谱峰更宽，分辨率更低；M=16时，相对于M=8的谱峰更窄，分辨率更高，但旁瓣数量也更多。（4）改变信噪比 SNR=0dB,20 dB ，其他输入同步骤 1，观察信噪比对谱估计性能的影响。SNR=0dB SNR=20 dB图5：周期图角度谱估计（信噪比）由图5可知，信噪比越大，所得角度谱谱峰比旁瓣要

8、越高，那么当信噪比小到一定程度，旁瓣会越来越高，容易与主瓣混淆，影响谱估计性能。3、编写Capon 最大似然法进行角度谱估计的程序。运行程序，输入参数同步骤1，观察谱峰位置是否正确，并重复步骤2 的内容。步骤1：如图6所示，谱峰位于-20°和30°处，与信号源的入射角度相符。图6：Capon角度谱估计步骤2：（1） d=4 d=2图7：Capon角度谱估计（阵元间距）对于空间采样，与奈奎斯特采样定理类似，必须满足d2，所以d=4时，谱峰位置正确，而d=2时，类似欠采样，除了在正确的角度上出现谱峰外，在其他多处地方也出现了谱峰，所以谱峰位置不正确。而由d=4与d=2时的谱峰对

9、比可知，d=2时分辨率更高，因此在满足d2的条件下，d越大，谱分辨率越高。（2）=0° =45°=90°图8：Capon角度谱估计（入射角度）由图8可知，在=0°,45°时，角度谱上在相应位置出现正确的谱峰，而在=90°时，在90°和-90°处都出现一个谱峰，影响DOA估计心能。在实际中，90°和-90°物理朝向相差明显，可以人工判断到底是90°还是-90°。（3）M=4 M=16图9：Capon角度谱估计（阵元数目）由图9可知，M=4时，相对于M=8的谱峰更宽，分辨率更低；

10、M=16时，相对于M=8的谱峰更窄，分辨率更高，Capon旁瓣数量比周期图旁瓣数量要少得多。（4）SNR=0dB SNR=20 dB图10：Capon角度谱估计（信噪比）由图10可知，信噪比越大，所得角度谱谱峰比旁瓣要越高，那么当信噪比小到一定程度，旁瓣会越来越高，容易与主瓣混淆，影响谱估计性能。4、编写MUSIC方法进行角度谱估计的程序。运行程序，输入参数同步骤1，观察谱峰位置是否正确，并重复步骤2 的内容。步骤1：如图11所示，谱峰位于-20°和30°处，与信号源的入射角度相符。图11：MUSIC角度谱估计步骤2：（1） d=4 d=2图12：MUSIC角度谱估计（阵元

11、间距）对于空间采样，与奈奎斯特采样定理类似，必须满足d2，所以d=4时，谱峰位置正确，而d=2时，类似欠采样，除了在正确的角度上出现谱峰外，在其他多处地方也出现了谱峰，所以谱峰位置不正确。而由d=4与d=2时的谱峰对比可知，d=2时分辨率更高，因此在满足d2的条件下，d越大，谱分辨率越高。（2） =0° =45°=90°图13：MUSIC角度谱估计（入射角度）由图13可知，在=0°,45°时，角度谱上在相应位置出现正确的谱峰，而在=90°时，在90°和-90°处都出现一个谱峰，影响DOA估计心能。（3） M=4 M

12、=16图14：MUSIC角度谱估计（阵元数目）由图14可知，M=4时，相对于M=8的谱峰更宽，分辨率更低；M=16时，相对于M=8的谱峰更窄，分辨率更高，MUSIC几乎没有旁瓣。（4）SNR=0dB SNR=20 dB图15：MUSIC角度谱估计（信噪比）由图15可知，信噪比越大，所得角度谱谱峰比旁瓣要越高，那么当信噪比小到一定程度，旁瓣会越来越高，容易与主瓣混淆，影响谱估计性能。5、改变两个源信号入射角度分别为-5°,5°，其他输入参数同步骤1，将周期图法、Capon 最大似然法以及MUSIC方法的角度谱估计结果画在同一坐标，比较三种方法在分辨率上的差别。如图16所示，入

13、射角度为-5°,5°，使用周期图法根本无法区分，分辨率太低。Capon最大似然法能够区分两个角度，但两峰有明显交界处，分辨率比周期图要高。MUSIC方法能够明显地区分两个角度，且交界处比较少，分辨率最高。所以三种方法分辨率比较为：MUSIC > Capon >周期图法。图16：周期图、Capon、MUSIC角度谱估计比较图四、 实验结论通过实验可得：使用周期图法、Capon法、MUSIC法都能在一定条件下实现DOA估计方法。阵元间距d的大小能够影响到角度谱分辨率，d越大，分辨率越高，但必须满足d2，否则由于采样率过低导致结果错误。 一般情况下，入射角对谱估计性能

14、没有影响，但当入射角接近90度或-90度，谱峰也会出现在-90度或90度处，影响判断。阵元数目越多，分辨率越高，对于周期图法，旁瓣越明显。信噪比越高，所得谱峰越明显，旁瓣越少。当信噪比低于一定程度，旁瓣会影响到谱估计结果。三种方法的分辨率不一样，其分辨率大小比较为：MUSIC > Capon >周期图法。五、 思考题1 推导周期图法进行角度谱估计的矩阵表达形式：PBT=IM=1MEHRE2 当入射源信号不是彼此独立的随机信号，而是频率相同的正弦信号（相干信源）时，上述三种角度谱估计方法的性能有什么变化？六、 附录DOA.mclose all;clc;clear;M = 8;lumb

15、t = 1;d = lumbt/2;N = 100;SNR =10;f = 300000000;th1 = -20/180*pi;th2 = 30/180*pi;c=3*108;% generate signalTable = (1/sqrt(2)*-1-1i -1+1i 1-1i 1+1i;T1 = randi(0,1,1,N*2);T_IQ1 = reshape(T1,2,N).'symbol1 = bi2de(T_IQ1,'left-msb');s1 = Table(symbol1+1);T2 = randi(0,1,1,N*2);T_IQ2 = reshape(

16、T2,2,N)'symbol2 = bi2de(T_IQ2,'left-msb');s2 = Table(symbol2+1);s = s1;s2;% generate white complex gaussian noisen = 1/(sqrt(2*10(SNR/10)*(randn(M,N) + 1i*randn(M,N);% Receive transformation matrixw = 2*pi*f;t1 = d*sin(th1)/c;a1 = exp(1i*w*(0:M-1)*t1);t2 = d*sin(th2)/c;a2 = exp(1i*w*(0:M

17、-1)*t2);A = a1;a2.'X = A*s + n;% PeriodogramR = X*X'/N;for k=0:1:180 th = pi/180*k-pi/2; tth = d*sin(th)/c; E = (exp(1i*w*tth*(0:M-1).' P_BT(k+1) = E'*R*E/M;endfigure(1),clfplot(-90:1:90,10*log10(abs(P_BT);% legend('Periodogram');%hold on% CaponR_inv = inv(R);for k=0:1:180 th

18、 = pi/180*k-pi/2; tth = d*sin(th)/c; E = (exp(1i*w*tth*(0:M-1).' P_Capon(k+1) = 1/(E'*R_inv*E);end figure(2),clfplot(-90:1:90,10*log10(abs(P_Capon);% legend('Capon');% MUSICMUSIC_N = n*n'/N;MUSIC_delt = mean(diag(MUSIC_N);MUSIC_S = s*s'/N;MUSIC_R = A*MUSIC_S*A' + MUSIC_N;V,Lbt = eig(MUSIC_R);Lbt1 = diag(Lbt);p = find(abs(Lbt1-MUSIC_delt)<0.1);for k=0:1:180 th =