基于MUSIC算法的测向性能仿真

基于MUSIC算法的测向性能仿真 ----- MUSIC 2013 年 1 月 16 日 ----- 摘 要 随着移动通信技术的飞速发展，智能天线技术研究的不断深入，来波方向(DOA)估计技术逐渐成为研究的热点之一，而MUSIC算法是智能天线技术的典型算法。本文在对MUSIC算法进行分析的基础上， 设计 领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计 了MUSIC算法的仿真程序，对不同情况下该算法的性能进行了仿真分析。仿真结果表明该算法在不同阵列结构、信号入射角度时具有不同的性能。 关键词:智能天线;DOA;MUSIC;阵元 ----- 目 录 摘 要................................................................................................................................................I 引 言............................................................................................................................................... 1 一、MUSIC算法介绍 ..................................................................................................................... 1 1.1 MUSIC算法的提出 ........................................................................................................... 1 1.2波达方向估计问题中的阵列信号数学模型 ..................................................................... 2 1.3阵列协方差矩阵的特征分解 ............................................................................................. 4 1.4 MUSIC算法的原理及实现 ............................................................................................... 6 1.5 MUSIC算法的实现步骤: ............................................................................................... 8 二、MUSIC算法的DOA估计仿真 .............................................................................................. 8 2.1MUSIC算法的基本仿真 .................................................................................................... 8 2.2 MUSIC算法DOA估计与阵元数的关系......................................................................... 9 2.3 MUSIC算法DOA估计与阵元间距的关系................................................................... 10 2.4 MUSIC算法DOA估计与快拍数的关系....................................................................... 11 2.5 MUSIC算法DOA估计与信噪比的关系....................................................................... 12 2.6 MUSIC算法DOA估计与信号入射角度差的关系 ....................................................... 13 三、MUSIC算法性能分析小结 ................................................................................................... 15 参考文献................................................................................................................................. 15 附 录..................................................................................................................................... 16 附录一:MUSIC算法的基本仿真源代码 ........................................................................... 16 附录二:MUSIC算法DOA估计与不同阵元数关系仿真源代码 ..................................... 17 附录三: MUSIC算法DOA估计与阵元间距的关系仿真源代码 ................................... 18 附录四:MUSIC算法DOA估计与快拍数的关系仿真源代码 ......................................... 21 附录五: MUSIC算法DOA估计与信噪比的关系仿真源代码 ....................................... 22 附录六:MUSIC算法DOA估计与信号入射角度差的关系仿真源代码 ......................... 24 图目录 图1-1 等距线阵与远场信号 .................................................................................................. 2 图2-1 MUSIC算法的DOA估计谱 ...................................................................................... 9 图2-2阵元数不同时MUSIC算法的DOA估计谱 ............................................................ 10 图2-3 阵元间距不同时MUSIC算法的DOA估计谱 ....................................................... 11 图2-4 快拍数不同时MUSIC算法的DOA估计谱 ........................................................... 12 图2-5 信噪比不同时MUSIC算法的DOA估计谱 ........................................................... 13 图2-6 角度间隔不同时MUSIC算法的DOA估计谱 ....................................................... 14 ----- 引 言 智能天线技术是当前无线移动通信领域颇为关注和研究的热点领域之一，可将无线电的信号导向到具体的方向上，产生空间定向波束，使天线主波束对准用户信号到达方向，旁瓣或零陷对准干扰信号的到达方向，起到充分高效利用移动用户信号并删除或抑制干扰信号的目的。而波束形成的关键是要准确知道信号的到达方向，即波达方向，所以波达角估计(DOA)是波束形成的基础。本文着重分析了用于DOA估计的典型算法--MUSIC(Multiple Signal Classification)算法，然后对不同的条件下MUSIC算法的性能进行了Matlab的仿真和分析。 一、MUSIC算法介绍 1.1 MUSIC算法的提出 多重信号分类(MUSIC)算法是Schmidt等人在1979年提出的。这一算法的提出开创了空间谱估计算法研究的新时代，促进了特征结构类算法的兴起和发展，该算法已成为空间谱估计理论体系中的标志性算法。此算法提出之前的有关算法都是针对阵列接收数据协方差矩阵进行直接处理，而MUSIC算法的基本思想则是对任意阵列输出数据的协方差矩阵进行特征分解，从而得到与信号分类相对应的信号子空间和与信号分量相正交的噪声子空间，然后利用这两个子空间的正交性构造空间谱函数，通过谱峰搜索，检测信号的DOA。 正是由于MUSIC算法在特定的条件下具有很高的分辨力、估计精度及稳定性，从而吸引了大量的学者对其进行深入的研究和分析。总的来说，它用于阵列的波达方向估计有以下一些突出的优点: (1)多信号同时测向能力 (2)高精度测向 (3)对天线波束内的信号的高分辨测向 (4)可适用于短数据情况 ----- (5)采用高速处理技术后可实现实时处理 1.2波达方向估计问题中的阵列信号数学模型 为了分析推导的方便，现将波达方向估计问题中的数学模型作理想状态的假设如下: (1)各待测信号源具有相同的极化、且互不相关的。一般考虑信号源为窄带的， ,且各信号源具有相同的中心频率。待测信号源的个数为D。 0 (2)天线阵列是由M(M>D)个阵元组成的等间距直线阵，各阵元特性相同，各向同性，阵元间隔为d，并且阵元间隔不大于最高频率信号半波长。 (3)天线阵列处于各信号源的远场中，即天线阵列接收从各信号源传来的信号为平面波。 2,(4)各阵元上有互不相关，与各待测信号也不相关，方差为的零均值高斯白噪n(t)声。 m (5)各接收支路具有完全相同的特性。 , k d 1 2 3 M 图1-1 等距线阵与远场信号 S(t)设由第k(k=1，2，…D)个信号源辐射到天线阵列的波前信号为，前k ----- S(t)S(t)面已假设为窄带信号，则可以表示为以下形式: kk S(t),s(t)exp{j,t} (1.1) kkk s(t)S(t),S(t)式中是的复包络，是信号的角频率。前面已经假设D个信号kkkk 具有相同的中心频率，所以有: ,2c,, (1.2) ,,k0, 式中c是电磁波波速，,是公用的信号波长。 t设电磁波通过天线阵列尺寸所需的时间为，则根据窄带假设，有如下近似: 1 S(t,t),S(t) (1.3) k1k 故延迟后的波前信号为: ~S(t,t),S(t,t)exp[j,(t,t)],S(t)exp[j,(t,t)]k1k101k01 (1.4) 所以，若以第一个阵元为参考点，则t时刻等间距直线阵中的第m(m=1，2，…M) 个阵元对第k个信号源的感应信号为: ,,2sindk()exp[(1)]aSt,jm,kk, (1.5) a其中，为第m个阵元对第k个信号源的影响，前面以假设各阵元无方向性，k dsin,ka,1,所以可取。为第k个信号源的方位角，m表示由第m个阵(,1)kkc元与第1个阵元间的波程差所引起的信号相位差。 计及测量噪声和所有信号源来波，第m个阵元的输出信号为: D,,2dsink (1.6)x(t),S(t)exp[,j(m,1)]，n(t),mkm,,1k n(t) 其中是测量噪声，所有标号为m表示该量属于第m个阵元，所有标号 m 为k表示该量属于第k个信号源。 设 ,,2sindk, (1.7) ()exp[(1)]a,,jm,mk, ----- 为第m个阵元对第k个信号源的响应函数。 则第m个阵元的输出信号为: D (1.8) x(t),a(,)S(t)，n(t),mmkkm,1k S(t)其中是第k个信号源在阵元上的信号强度。 k 运用矩阵的定义，可以得到更为简洁的表达式: X=AS+N (1.9) 式中 T 1.10 ()X,[x(t),x(t),...x(t)]12M TS,[S(t),S(t),...S(t)] 1.11 ()12D T A,[a(,),a(,),...a(,)]12D = 1.12 ()11...1，，,,,,,j,j,j,12Dee...e,, ,,............ ,,,j(M,1,),j(M,1,),j(M,1),12Dee...e，， 1.13 (),2d,,sin,kk, T 1.14 ()N,[n(t),n(t),...,n(t)]12M x(t)x(t)对进行N点采样，要处理的问题就变成了通过输出信号的采样mm{x(i),i,1,2,...,M},,,...,,估计出信号源的波达方向角。 m12,D 由此，可以很自然的将阵列信号看作是噪声干扰的若干空间谐波的叠加，从 而将波达方向估计问题与谱估计联系起来。 1.3阵列协方差矩阵的特征分解 R对阵列输出x作相关处理，得到其协方差矩阵: x HR,E[XX] (1.15) x ----- 其中，H表示矩阵共轭转置。 前面已假设信号与噪声互不相关、且噪声为零均值白噪声，因此将式(1.9)代入式(1.15)，可以得到: HR,E[(AS，N)(AS，N)]x HHHAE[SS]A，E[NN]= HARA，R = (1.16) sN 式中 HR,E[SS] (1.17) s 称为信号的相关矩阵。 2R,,I (1.18) N 2,是噪声的相关矩阵，是噪声功率，I是M*M阶的单位矩阵。 :RR实际应用中，通常无法直接得到，能使用的只有样本的协方差矩阵: xx N:1H (1.19) R,X(i)X(i),xN,1i :RRN,,是的最大似然估计，当采样数时，它们是一致的，但实际情况xx 中将由于样本数有限而造成误差。 根据矩阵特征分解的理论，可以对阵列协方差矩阵进行特征分解。首先考虑理想情况，即无噪声的情况: HR,ARA (1.20) xs 对于均匀线阵，矩阵A是由式(1.12)所定义的范德蒙德矩阵，只要满足: ,,,i,j (1.21) ij RRank(R),D则，它的各列相互独立，这样，若为非奇异矩阵(，各信号源ss两两不相干)，且M>D，则有: Hrank(ARA),D (1.22) s HR,E[XX]由于，所以有: x ----- HR,R (1.23) xx HRRARA即是Hermite矩阵，它的特征值都是实数。又由于是正定的，因此矩阵xss是半正定的，它有D个正特征值和M-D个零特征值。 再考虑有噪声存在的情况 H2R,ARA，,I (1.24) xs 2RR,,,,,...,,由于>0，为满秩阵，所以有M个正实特征值，分别对xx12M Rv,v,...v应于M个特征向量。又由于是Hermite矩阵，所以各特征向量是相x12M 互正交的，即: Hvv,0i,j (1.25) ij H2ARA,与信号有关的特征值只有D个，分别等于矩阵的各特征值与之s 22,,和，其余的M-D个特征值为，也就是说，是R的最小特征值，它是M-D v维的。对应的特征向量，i=1,2，…，M中，也有D个是与信号有关的，另外i M-D个是与噪声有关的，在下一节里，将利用以上这些特征分解的性质求出信 ,号源的波达方向。 k 1.4 MUSIC算法的原理及实现 通过对阵列协方差矩阵的特征分解，可以得到如下结论: R将矩阵的特征值进行从小到大的排序，即 x ,,,,...,,,0 (1.26) 12M 其中D个较大的特征值对应于信号，M-D个较小的特征值对应于噪声。 R矩阵的属于这些特征值的特征向量也分别对应于信号和噪声，因此，可x R以把的特征值(特征向量)划分为信号特征值(特征向量)与噪声特征值(特征向x 量)。 ,Rv,设是矩阵的第i个特征值，是与个相对应的特征向量，则有: ixii ----- Rv,,v (1.27) xiii 2,,,R再设是的最小特征值 ix 2Rv,,v i=D+1，D+2，…M (1.28) xii 将 H2R,ARA，,I (1.29) xs 代入上式，可得: H22,v,(ARA，,I)v (1.30) isi将上式右边展开与左边比较，可得: HARAv,0 (1.31) si ,1HH,1AAR(AA)因是D*D维的满秩矩阵，存在;而同样存在，则上式两边同s ,1H,1HR(AA)A乘以后变成: s ,1H,1HHR(AA)AARAv,0 (1.32) ssi于是有 HAv,0 i=D+1，D+2，…，M (1.33) i v上式表明:噪声特征值所对应的特征向量(称噪声特征向量)，与矩阵A的i列向量正交，而A的各列是与信号源的方向相对应的。这就是利用噪声特征向 量求解信号源方向的出发点。 E用各噪声特征向量为列，构造一个噪声矩阵: n E,[v,v,...,v] (1.34) nD，1D，2M P(,)定义空间谱: mu 11, P() = (1.35) ,mu2HHHa(,)EEa(,)Ea(,)nnn Ea(,)该式中分母是信号向量和噪声矩阵的内积，当和的各列正交时，该n P(,)分母为零，但由于噪声的存在，它实际上为一最小值，因此有一尖峰。由mu ----- 该式，使变化，通过寻找波峰来估计到达角。 , 1.5 MUSIC算法的实现步骤: (1)根据N个接收信号矢量得到下面协方差矩阵的估计值: N1H (1.36) R,X(i)X(i),xN,1i 对上面得到的协方差矩阵进行特征值分解 H2R,ARA，,I (1.37) xs (2)按特征值的大小顺序，把与信号个数D相等的特征值和对应的特征向量看作信号部分空间，把剩下的M-D个特征值和特征向量看作噪声部分空间。得 E到噪声矩阵: n HAv,0 i=D+1，D+2，…，M (1.38) i E,[v,v,...,v] (1.39) nD，1D，2M ,(3)使变化，按照式 1,P() (1.40) ,muHHa(,)EEa(,)nn 来计算谱函数，通过寻求峰值来得到波达方向的估计值。 二、MUSIC算法的DOA估计仿真 2.1MUSIC算法的基本仿真 模拟2个独立窄带信号分别以20?，70?的方向入射到均匀线阵上，信号间互不相关，与噪声相互独立，噪声为理想高斯白噪声，阵元间距为入射信号波长的1,2，快拍数为200，信噪比为20dB，阵元数分别为10，50，100。其仿真结果如图2-1所示: ----- MUSIC算法的DOA估计谱 0 -10 -20 ) /dB, -30 谱函数P(-40 -50 -60-100-80-60-40-20020406080100 ,/degree角度 图2-1 MUSIC算法的DOA估计谱 由图2-1可以看出在符合假设的前提下，采用MUSIC算法能构造出针状的谱峰，可以很好的估计出入射信号的个数和方向，能有效的估计出独立信号源的DOA,并且在模型准确的前提下，对DOA的估计可以达到任意精度，克服了传统测向定位方法精度低的缺点 ,可以有效解决密集信号环境中多个辐射源的高分辨率、高精度测向定位问题。可以看出超分辨率的 MUSIC算法具有测向准确度、灵敏度高的特点且具有潜在分辨多信号的能力，具有较好的性能和较高的效率，能提供高分辨率及渐近无偏的到达角估计，这对实际中的应用具有十分重要的意义。 2.2 MUSIC算法DOA估计与阵元数的关系 模拟2个独立窄带信号分别以20?，60?的方向入射到均匀线阵上，信号间互不相关，与噪声相互独立，噪声为理想高斯白噪声，阵元间距为入射信号波长的1,2，快拍数为200，信噪比为20dB，阵元数分别为5，10，20。其仿真结果如图2-2所示: ----- MUSIC算法的DOA估计谱 0 阵元数为5 阵元数为10 阵元数为20-10 -20 ) /dB, -30 谱函数P(-40 -50 -60 -100-80-60-40-20020406080100 ,/degree角度 图2-2阵元数不同时MUSIC算法的DOA估计谱 由图2-2可以看出，其他条件不变的情况下，随着阵元数的增加，DOA估计谱的波束宽度变窄，阵列的指向性变好，也就是说阵列分辨空间信号的能力增强。由此可以看出，要得到更加精确的DOA估计谱，可以增加阵元数量，但阵元数量越多，需要处理的数据越多，运算量越大，运行速度越慢。由上图可以看出阵元数大到一定数量时，波形变化不会很明显。因此，在实际应用中，可根据具体条件适当选取阵元数量，在确保估计谱准确的前提下，尽量减少资源浪费，加快运行的速度，提高工作效率。 2.3 MUSIC算法DOA估计与阵元间距的关系 模拟2个独立窄带信号分别以20?，60?的方向入射到均匀线阵上，信号间互不相关，与噪声相互独立，噪声为理想高斯白噪声，阵元数为10，快拍数 ,,,为200，信噪比为20dB，阵元间距分别为/6、/2、。其仿真结果如图2-3所示: ----- MUSIC算法的DOA估计谱 0 间距为六分之波长 间距为半波长 间距为一个波长-10 -20 ) /dB, -30 谱函数P(-40 -50 -60 -100-80-60-40-20020406080100 ,/degree角度 图2-3 阵元间距不同时MUSIC算法的DOA估计谱 由图2-3可以看出，在其他条件不变的前提下，当阵元间距不大于半波长时，随着阵元间距的增加，DOA估计谱的波束宽度变窄，阵列的指向性变好，也就是说MUSIC算法的分辨力随着阵元间距的加大相应提高，但当阵元间距大于半波长时，估计谱除了信号源方向外在其他方向出现了虚假谱峰，也就失去了估计的准确性。可见，在实际应用中，要十分注意阵元间的距离，可以适当增加阵元间距但绝不能超过半波长，这一点非常重要，最好是将阵元间距设为半波长。 2.4 MUSIC算法DOA估计与快拍数的关系 模拟2个独立窄带信号分别以20?，60?的方向入射到均匀线阵上信号间互不相关，与噪声相互独立，噪声为理想高斯白噪声，阵元数为10，阵元间距为入射信号波长的1,2，信噪比为20dB，快拍数分别为5，50，200。其仿真结果如图2-4所示: ----- MUSIC算法的DOA估计谱 0 快拍数为5 快拍数为50 快拍数为200-10 -20 ) /dB, -30 谱函数P(-40 -50 -60 -100-80-60-40-20020406080100 ,/degree角度 图2-4 快拍数不同时MUSIC算法的DOA估计谱 由图2-4可以看出，在其他条件不变的情况下，随着快拍数的增加，DOA估计谱的波束宽度变窄，阵列的指向性变好，阵列分辨空间信号的能力增强，MUSIC算法的估计精度增加。由此可见，可通过增加采样快拍数来增加DOA估计的精确度，但是采样快拍数越多，需要处理的数据就越多，MUSIC算法的运算量就越大，速度就越慢，所以在实际应用中要合理的选取采样快拍数，在确定DOA估计谱准确的前提下，尽量减少运算量，加快工作速度，节省人力物力，节约资源。 2.5 MUSIC算法DOA估计与信噪比的关系 模拟2个独立窄带信号分别以20?，60?的方向入射到均匀线阵上，信号间互不相关，与噪声相互独立，噪声为理想高斯白噪声，阵元数为10，阵元间距为入射信号波长的1,2，快拍数为200，信噪比分别为-10dB，0dB，20dB。其仿真结果如图2-5所示: ----- MUSIC算法的DOA估计谱 0 信噪比为-10dB 信噪比为0dB 信噪比为20dB-10 -20 ) /dB, -30 谱函数P(-40 -50 -60 -100-80-60-40-20020406080100 ,/degree角度 图2-5 信噪比不同时MUSIC算法的DOA估计谱 由图2-5可以看出，在其他条件不变的情况下，随着信噪比的增加，DOA估计谱的波束宽度变窄，阵列的指向性变好，MUSIC算法的分辨力增加，信噪比的高低直接影响着超分辨方位估计算法的性能。在低信噪比时，MUSIC算法的性能会急剧下降，因而提高算法在低信噪比条件下的估计性能是超分辨DOA算法的研究重点。有学者提出了一种基于多级维纳滤波器(MSWF:Multistage Weiner Filtering)的信号波达方向(DOA)估计算法，该算法在信号可能入射方向用MSWF 估计信号子空间，并在 MSWF 分解后互相关函数最小以及信号子空间估值与噪声子空间正交时判定估计有效，进而构造空间谱来实现信号 DOA 估计。已证明在低信噪比条件下，该算法比子空间类算法有更好的分辨率和误差性能。低信噪比条件下对DOA的精确估计还有很大的发展改进空间，有待进一步研究。 2.6 MUSIC算法DOA估计与信号入射角度差的关系 模拟2个独立窄带信号入射到均匀线阵上，信号间互不相关，与噪声相互独 ----- 立，噪声为理想高斯白噪声，阵列的阵元数为10，快拍数为200。阵元间距为入射信号波长的1,2，信噪比为20dB，信号入射角度差分别为5?，10?，40?。其仿真结果如图2-6所示: MUSIC算法的DOA估计谱 0 角度差为5? 角度差为10? 角度差为40?-10 -20 ) /dB, -30 谱函数P(-40 -50 -60 -100-80-60-40-20020406080100 ,/degree角度 图2-6 角度间隔不同时MUSIC算法的DOA估计谱 图2-6说明在其他条件不变的情况下，随着信号入射角度差的增加，DOA估计谱的波束宽度变窄，阵列的指向性变好，MUSIC算法的分辨力增加。当信号来波方向间隔角度很小时，不能准确估计信号源数。通常的阵列信号源数估计方法, 都是在信号来波方向角度差较大情况下进行的估计，当信号来波方向的角度差比较小时, 这些方法估计都要失效。已有学者提出了平方根修正的 Gerschgorin半径估计方法，对信号来波方向角度差小时,也能很好的估计信号源数。目前提出的一些信号源数估计方法，大都存在一定的应用条件，这在一定程度上限制了DOA算法的实际运用，因此，研究符合实际应用环境的实时、稳健的信号源数与DOA的联合估计等仍具有十分重要的现实意义。 ----- 三、MUSIC算法性能分析小结 通过上述几组仿真可以看出超分辨率的 MUSIC算法具有较好的性能和较高的效率，能提供高分辨率及渐近无偏的到达角估计。而且阵元数越多，快拍数越多，信噪比越高，信号入射角度差越大 MUSIC算法的分辨率越高，当阵元间距不大于载波半波长时， MUSIC算法的分辨力随着阵元间距的加大相应提高， /2时，空间谱除了信号源方向外在其他方向出现虚假谱峰。但当阵元间距大于, 在小信噪比和小的角度间隔时，MUSIC算法的估计性能下降，已有学者提出了一些的改进算法，但这些问题仍是当前研究的热点。此外，有学者提出了很多对MUSIC算法的改进。可见，MUSIC算法还有很大的发展空间，值得我们进一步研究。 参考文献 [1]Eli张贤达，保铮(通信信号处理[M](北京:国防工业出版社，2000( [2]张贤达(现代信号处理[M](北京:清华大学出版社，2002( [3]郑洪(MUSIC算法与波达方向估计研究[D](成都:四川大学，2005( [4]张娟(智能天线中DOA估计算法的研究[D](哈尔滨:哈尔滨工业大学，2006( [5]张玲华，郑宝玉(随机信号处理[M](北京:清华大学出版社，2003( [6]Simon Haykin(自适应滤波器原理[M](北京:电子工业出版社，2003( [7]楼顺天，李博菡. 基于MATLAB的系统分析与设计——信号处理[M]. 西安: 西安电子科技大学出版社, 2000. [8]王正林，刘明.精通MATLAB7[M].北京:电子工业出版社，2006. [9]王进，赵拥军，王志刚.低信噪比条件下的高分辨DOA估计算法[J].计算机工程，2005.Vol.35.No.3. ----- 附 录 附录一:MUSIC算法的基本仿真源代码 clc clear all format long %将数据显示为长整型科学计数 N=200;%快拍数 doa=[20 60]/180*pi; %信号到达角 w=[pi/4 pi/3]';%信号频率 M=10;%阵元数 P=length(w); %信号个数 lambda=150;%波长 d=lambda/2;%阵元间距 snr=20;%信噪比 B=zeros(P,M); %创建一个P行M列的0矩阵 for k=1:P B(k,:)=exp(-j*2*pi*d*sin(doa(k))/lambda*[0:M-1]); %矩阵赋值 end B=B'; xx=2*exp(j*(w*[1:N])); %仿真信号 x=B*xx; x=x+awgn(x,snr);%加入高斯白噪声 R=x*x'; %数据协方差矩阵 [U,V]=eig(R); %求R的特征值和特征向量 UU=U(:,1:M-P); %估计噪声子空间 theta=-90:0.5:90; %%谱峰搜索 for ii=1:length(theta) AA=zeros(1,length(M)); for jj=0:M-1 AA(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda); end WW=AA*UU*UU'*AA'; Pmusic(ii)=abs(1/ WW); end Pmusic=10*log10(Pmusic/max(Pmusic)); %空间谱函数 plot(theta,Pmusic,'-k') xlabel('角度 \theta/degree') ylabel('谱函数P(\theta) /dB') title('MUSIC算法的DOA估计谱') grid on ----- 附录二:MUSIC算法DOA估计与不同阵元数关系仿真源代码 clc clear all format long %将数据显示为长整型科学计数 N=200; %快拍数 doa=[20 60]/180*pi; %信号到达角 w=[pi/4 pi/3]'; %信号频率 M1=5; %阵元数 M2=10; M3=20; P=length(w); %信号个数 lambda=150; %波长 d=lambda/2; %阵元间距 snr=20; %信噪比 B1=zeros(P,M1); B2=zeros(P,M2); B3=zeros(P,M3); for k=1:P B1(k,:)=exp(-j*2*pi*d*sin(doa(k))/lambda*[0:M1-1]); %矩阵赋值 B2(k,:)=exp(-j*2*pi*d*sin(doa(k))/lambda*[0:M2-1]); B3(k,:)=exp(-j*2*pi*d*sin(doa(k))/lambda*[0:M3-1]); end B1=B1'; B2=B2'; B3=B3'; xx=2*exp(j*(w*[1:N])); %仿真信号 x1=B1*xx; x2=B2*xx; x3=B3*xx; x1=x1+awgn(x1,snr); %加入高斯白噪声 x2=x2+awgn(x2,snr); x3=x3+awgn(x3,snr); R1=x1*x1'; %数据协方差矩阵 R2=x2*x2'; R3=x3*x3'; [U1,V1]=eig(R1); %求R的特征值和特征向量 [U2,V2]=eig(R2); [U3,V3]=eig(R3); UU1=U1(:,1:M1-P); ; %估计噪声子空间 UU2=U2(:,1:M2-P); UU3=U3(:,1:M3-P); theta=-90:0.5:90; ----- %%谱峰搜索 for ii=1:length(theta) AA1=zeros(1,length(M1)); for jj=0:M1-1 AA1(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda); end WW1=AA1*UU1*UU1'*AA1'; Pmusic1(ii)=abs(1/ WW1); end for ii=1:length(theta) AA2=zeros(1,length(M2)); for jj=0:M2-1 AA2(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda); end WW2=AA2*UU2*UU2'*AA2'; Pmusic2(ii)=abs(1/ WW2); end for ii=1:length(theta) AA3=zeros(1,length(M3)); for jj=0:M3-1 AA3(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda); end WW3=AA3*UU3*UU3'*AA3'; Pmusic3(ii)=abs(1/ WW3); end Pmusic1=10*log10(Pmusic1/max(Pmusic1)); %空间谱函数 Pmusic2=10*log10(Pmusic2/max(Pmusic2)); Pmusic3=10*log10(Pmusic3/max(Pmusic3)); plot(theta,Pmusic1,'--k','LineWidth',2.0) hold on plot(theta,Pmusic2,'r','LineWidth',2.0) hold on plot(theta,Pmusic3,'-b','LineWidth',2.0) hold off xlabel('角度 \theta/degree') ylabel('谱函数P(\theta) /dB') title('MUSIC算法的DOA估计谱') grid on 附录三: MUSIC算法DOA估计与阵元间距的关系仿真源代码 clc clear all format long %将数据显示为长整型科学计数 ----- N=200;%快拍数 doa=[20 60]/180*pi; %信号到达角 w=[pi/4 pi/3]';%信号频率 M=10;%阵元数 P=length(w); %信号个数 lambda=150;%波长 d=lambda/6;%阵元间距 snr=20;%信噪比 B=zeros(P,M); %创建一个P行M列的0矩阵 for k=1:P B(k,:)=exp(-j*2*pi*d*sin(doa(k))/lambda*[0:M-1]); %矩阵赋值 end B=B'; xx=2*exp(j*(w*[1:N])); %仿真信号 x=B*xx; x=x+awgn(x,snr);%加入高斯白噪声 R=x*x'; %数据协方差矩阵 [U,V]=eig(R); %求R的特征值和特征向量 UU=U(:,1:M-P); %估计噪声子空间 theta=-90:0.5:90; %%谱峰搜索 for ii=1:length(theta) AA=zeros(1,length(M)); for jj=0:M-1 AA(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda); end WW=AA*UU*UU'*AA'; Pmusic(ii)=abs(1/ WW); end Pmusic=10*log10(Pmusic/max(Pmusic)); %空间谱函数 plot(theta,Pmusic ,'--k','linewidth',2.0) hold on clc clear all format long %将数据显示为长整型科学计数 N=200;%快拍数 doa=[20 60]/180*pi; %信号到达角 w=[pi/4 pi/3]';%信号频率 M=10;%阵元数 P=length(w); %信号个数 lambda=150;%波长 d=lambda/2;%阵元间距 snr=20;%信噪比 B=zeros(P,M); %创建一个P行M列的0矩阵 ----- for k=1:P B(k,:)=exp(-j*2*pi*d*sin(doa(k))/lambda*[0:M-1]); %矩阵赋值 end B=B'; xx=2*exp(j*(w*[1:N])); %仿真信号 x=B*xx; x=x+awgn(x,snr);%加入高斯白噪声 R=x*x'; %数据协方差矩阵 [U,V]=eig(R); %求R的特征值和特征向量 UU=U(:,1:M-P); %估计噪声子空间 theta=-90:0.5:90; %%谱峰搜索 for ii=1:length(theta) AA=zeros(1,length(M)); for jj=0:M-1 AA(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda); end WW=AA*UU*UU'*AA'; Pmusic(ii)=abs(1/ WW); end Pmusic=10*log10(Pmusic/max(Pmusic)); %空间谱函数 plot(theta,Pmusic,'r','linewidth',2.0) hold on clc clear all format long %将数据显示为长整型科学计数 N=200;%快拍数 doa=[20 60]/180*pi; %信号到达角 w=[pi/4 pi/3]';%信号频率 M=10;%阵元数 P=length(w); %信号个数 lambda=150;%波长 d=lambda;%阵元间距 snr=20;%信噪比 B=zeros(P,M); %创建一个P行M列的0矩阵 for k=1:P B(k,:)=exp(-j*2*pi*d*sin(doa(k))/lambda*[0:M-1]); %矩阵赋值 end B=B'; xx=2*exp(j*(w*[1:N])); %仿真信号 x=B*xx; x=x+awgn(x,snr);%加入高斯白噪声 R=x*x'; %数据协方差矩阵 [U,V]=eig(R); %求R的特征值和特征向量 ----- UU=U(:,1:M-P); %估计噪声子空间 theta=-90:0.5:90; %%谱峰搜索 for ii=1:length(theta) AA=zeros(1,length(M)); for jj=0:M-1 AA(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda); end WW=AA*UU*UU'*AA'; Pmusic(ii)=abs(1/ WW); end Pmusic=10*log10(Pmusic/max(Pmusic)); %空间谱函数 plot(theta,Pmusic,'-.b','linewidth',2.0) hold off xlabel('角度 \theta/degree') ylabel('谱函数P(\theta) /dB') title('MUSIC算法的DOA估计谱') grid on 附录四:MUSIC算法DOA估计与快拍数的关系仿真源代码 clc clear all format long %将数据显示为长整型科学计数 N1=5;%快拍数 N2=50; N3=200; doa=[20 60]/180*pi; %信号到达角 w=[pi/4 pi/3]';%信号频率 M=10;%阵元数 P=length(w); %信号个数 lambda=150;%波长 d=lambda/2;%阵元间距 snr=20;%信噪比 B=zeros(P,M); %创建一个P行M列的0矩阵 for k=1:P B(k,:)=exp(-j*2*pi*d*sin(doa(k))/lambda*[0:M-1]); %矩阵赋值 end B=B'; xx1=2*exp(j*(w*[1:N1])); %仿真信号 xx2=2*exp(j*(w*[1:N2])); xx3=2*exp(j*(w*[1:N3])); x1=B*xx1; x2=B*xx2; ----- x3=B*xx3; x1=x1+awgn(x1,snr);%加入高斯白噪声 x2=x2+awgn(x2,snr); x3=x3+awgn(x3,snr); R1=x1*x1'; %数据协方差矩阵 R2=x2*x2'; R3=x3*x3'; [U1,V1]=eig(R1); [U2,V2]=eig(R2); [U3,V3]=eig(R3); UU1=U1(:,1:M-P); UU2=U2(:,1:M-P); UU3=U3(:,1:M-P); theta=-90:0.5:90; %%谱峰搜索 for ii=1:length(theta) AA=zeros(1,length(M)); for jj=0:M-1 AA(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda); end WW1=AA*UU1*UU1'*AA'; WW2=AA*UU2*UU2'*AA'; WW3=AA*UU3*UU3'*AA'; Pmusic1(ii)=abs(1/ WW1); Pmusic2(ii)=abs(1/ WW2); Pmusic3(ii)=abs(1/ WW3); end Pmusic1=10*log10(Pmusic1/max(Pmusic1)); Pmusic2=10*log10(Pmusic2/max(Pmusic2)); Pmusic3=10*log10(Pmusic3/max(Pmusic3)); plot(theta,Pmusic1,'--k','linewidth',2.0) hold on plot(theta,Pmusic2,'-r','linewidth',2.0) hold on plot(theta,Pmusic3,'-.b','linewidth',2.0) hold off xlabel('角度 \theta/degree') ylabel('谱函数P(\theta) /dB') title('MUSIC算法的DOA估计谱') grid on 附录五: MUSIC算法DOA估计与信噪比的关系仿真源代码 clc ----- clear all format long %将数据显示为长整型科学计数 N=200;%快拍数 doa=[20 60]/180*pi; %信号到达角 w=[pi/4 pi/3]';%信号频率 M=10;%阵元数 P=length(w); %信号个数 lambda=150;%波长 d=lambda/2;%阵元间距 snr1=-10;%信噪比 snr2=0; snr3=20; B=zeros(P,M); %创建一个P行M列的0矩阵 for k=1:P B(k,:)=exp(-j*2*pi*d*sin(doa(k))/lambda*[0:M-1]); %矩阵赋值 end B=B'; xx=2*exp(j*(w*[1:N])); %仿真信号 x=B*xx; x1=x+awgn(x,snr1);%加入高斯白噪声 x2=x+awgn(x,snr2); x3=x+awgn(x,snr3); R1=x1*x1'; %数据协方差矩阵 R2=x2*x2'; R3=x3*x3'; [U1,V1]=eig(R1); [U2,V2]=eig(R2); [U3,V3]=eig(R3); UU1=U1(:,1:M-P); UU2=U2(:,1:M-P); UU3=U3(:,1:M-P); theta=-90:0.5:90; %%谱峰搜索 for ii=1:length(theta) AA=zeros(1,length(M)); for jj=0:M-1 AA(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda); end WW1=AA*UU1*UU1'*AA'; WW2=AA*UU2*UU2'*AA'; WW3=AA*UU3*UU3'*AA'; Pmusic1(ii)=abs(1/ WW1); Pmusic2(ii)=abs(1/ WW2); Pmusic3(ii)=abs(1/ WW3); ----- end Pmusic1=10*log10(Pmusic1/max(Pmusic1)); Pmusic2=10*log10(Pmusic2/max(Pmusic2)); Pmusic3=10*log10(Pmusic3/max(Pmusic3)); plot(theta,Pmusic1,'--k','linewidth',2.0) hold on plot(theta,Pmusic2,'-r','linewidth',2.0) hold on plot(theta,Pmusic3,'-.b','linewidth',2.0) hold off xlabel('角度 \theta/degree') ylabel('谱函数P(\theta) /dB') title('MUSIC算法的DOA估计谱') grid on 附录六:MUSIC算法DOA估计与信号入射角度差的关系仿真源代码 clc clear all format long %将数据显示为长整型科学计数 N=200;%快拍数 doa1=[20 25]/180*pi; %信号到达角 doa2=[20 30]/180*pi; doa3=[20 60]/180*pi; w=[pi/4 pi/3]';%信号频率 M=10;%阵元数 P=length(w); %信号个数 lambda=150;%波长 d=lambda/2;%阵元间距 snr=20;%信噪比 B=zeros(P,M); %创建一个P行M列的0矩阵 for k=1:P B1(k,:)=exp(-j*2*pi*d*sin(doa1(k))/lambda*[0:M-1]); %矩阵赋值 B2(k,:)=exp(-j*2*pi*d*sin(doa2(k))/lambda*[0:M-1]); B3(k,:)=exp(-j*2*pi*d*sin(doa3(k))/lambda*[0:M-1]); end B1=B1'; B2=B2'; B3=B3'; xx=2*exp(j*(w*[1:N])); %仿真信号 x1=B1*xx; x2=B2*xx; x3=B3*xx; ----- x1=x1+awgn(x1,snr);%加入高斯白噪声 x2=x2+awgn(x2,snr); x3=x3+awgn(x3,snr); R1=x1*x1'; %数据协方差矩阵 R2=x2*x2'; R3=x3*x3'; [U1,V1]=eig(R1); [U2,V2]=eig(R2); [U3,V3]=eig(R3); UU1=U1(:,1:M-P); UU2=U2(:,1:M-P); UU3=U3(:,1:M-P); theta=-90:0.5:90; %%谱峰搜索 for ii=1:length(theta) AA=zeros(1,length(M)); for jj=0:M-1 AA(1+jj)=exp(-j*2*jj*pi*d*sin(theta(ii)/180*pi)/lambda); end WW1=AA*UU1*UU1'*AA'; WW2=AA*UU2*UU2'*AA'; WW3=AA*UU3*UU3'*AA'; Pmusic1(ii)=abs(1/ WW1); Pmusic2(ii)=abs(1/ WW2); Pmusic3(ii)=abs(1/ WW3); end Pmusic1=10*log10(Pmusic1/max(Pmusic1)); Pmusic2=10*log10(Pmusic2/max(Pmusic2)); Pmusic3=10*log10(Pmusic3/max(Pmusic3)); plot(theta,Pmusic1,'--k','linewidth',2.0) hold on plot(theta,Pmusic2,'-r','linewidth',2.0) hold on plot(theta,Pmusic3,'-.b','linewidth',2.0) hold off xlabel('角度 \theta/degree') ylabel('谱函数P(\theta) /dB') title('MUSIC算法的DOA估计谱') grid on -----