阵列式信号处理MUSIC算法仿真yes

1、阵列式信号处理 MUSIC算法仿真渤海船舶职业学院常亮摘要：波达方向(DOA)f古计是阵列信号处理研究的重要内容之一多重信号分类 (MUSIC算法是一种经典有效的DOA方法。概述了阵列信号处理研究的主要内 容，详细介绍了多重信号分类算法的原理，并用 MATLAB进行仿真实现该算法， 研究了 MUSIC算法的性质。关键词：阵列信号处理多重信号分类算法仿真1、实验目的(1)掌握MUSIC算法进行DOA估计的基本原理(2) 了解影响DOA估计性能的因素2、实验内容假设阵元个数M=16,信号方向分别为0度、3度和-35度，噪声为高斯噪声 仿真：(1)快拍次数一定，分辨力与 SNR的关系(2)快拍次数一

2、定，分辨力与 孔径的关系；(3) SNR一定，分辨力与快拍的关系。3、实验原理窄带远场信号的DOA数学模型为X (t尸A(6 )S(t-N (t )阵列数据的协方差矩阵为R = E|XX H =A E |SS H A H + ；，I=ARs A H +；，I由于信号与噪声相互独立，数据协方差矩阵可以分解为与信号、噪声相关的两部分，其中Rs是信号的协方差矩阵，ARs AH是信号部分。对R进行特征分解有R=U SESU H +U nEn U H S S S N N N式中，U S是由大特征值对应的特征矢量张成的子空间也即信号子空间,UN而是由小特征值对应的特征矢量张成的子空间也即噪声子空间。根据前

3、面所述的性质2可知，在理想的条件下数据空间中的信号子空间与噪 声子空间是相互正交的，即信号空间中的导向矢量也与噪声子空间正交aH r UN =0经典的MUSIC算法正是基于上述这个性质提出的，但考虑到实际接收数据矩 阵是有限快拍数的，即数据协方差矩阵的最大似然估计为R = L xx HL i=1对!？进行特征分解可以计算得到噪声子空间特征矢量矩阵l?N 0由于噪声的存在，a(8)与巴并不能完全的正交。因此，实际上求 DOA是以最小优化搜索实现的，即 为usic =argmin a H F U N U N a所以，MUSIC算法的谱估计公式为Pmusic卜面给出MUSIC算法的计算步骤:(D 由

4、阵列的接收数据得到数据协方差矩阵 R ；对R进行特征分解；由R的特征值进行信号源数的判断；确定信号子空间Us与噪声子空间Un；根据信号参数的范围进行谱峰搜索；找出极大值点对应的角度就是信号入射方向。对于上述的MUSIC算法，还应该要注意以下几点：对于非理想情况下得到协方差矩阵的特征值满足下式:11 ,2 , -N , N 1 , - M所以判断信号源数需要用到有关信号源估计的方法来进行信号源数的确定。线阵的信号参数搜索范围为-90：90。，而面阵的搜索范围为-180：180口另外，一个确定阵列的导向矢量由阵元的位置唯一的确定；MUSIC算法的一种归一化形式，即caa FPMUSIC _ F -

5、 H -a 1 U N U N a 口在实际应用中，对于一维导向矢量有下式成立：aH 丁 a 二-M4、实验仿真及结果分析MUSIC随信噪比变化情况信号数共3个，分别位于2 , 30。和31。，它们功率相同，阵元数共16 个。这里快拍数取500,以保证快拍数足够，可以较好得得出统计特性，以去除 采样不足带来的分辨力低下。信噪比分别取为SNR=0,5,10,15,20dB得到图2。-10-20-30-40-50-1 0 0 - 80-60- 40-20020406080100从图中可以看到，随信噪比增大，MUSIC谱图旁瓣变低，主瓣变得更加尖锐，分 辨力变得更好。-5-10-15-20-25-3

6、011111-100-80-60-40-20020406080100信噪比为0dB信噪比为30dB TOC o 1-5 h z 1。I20,.30 I40,.50 , I60 ,.70 .80 I_.1111-100-80-60-40-20020406080100-1-2-3-4-5-6-770-8-100信噪比为50dB信噪比为-20dB信噪比为0dB时经典music算法解出的来波角度为 2.0700 , 31.7000 ,求 根music算法解出的来波角度为 30.5268 , 1.9996信噪比为30dB时经典music算法解出的来波角度为 2.0700 , 30.9900 , 32.3

7、200 ,求根 music算法解出的来波角度为 29.8490 , 1.9998 , 31.1206 。信噪比为50dB时经典 music算法解出的来波角度为2.0700 , 31 ,32.3100，求根 music算法解出的来波角度为 29.8602 ,2.0001 ,31.1104 。信噪比为-20dB时经典music算法和求根music算法都无法解出来波角度。由上面的图形可以看出，信噪比越大，越有利于信号角度的分辨，以及分辨 出角度差很小的两个信号，表现为峰值十分尖锐，而随着信噪比的降低，分辨力 下降，当为-20dB时已经无法分辨出30和31的两个信号，更糟糕的是还出 现了很多假的峰值，

8、导致对信号源数的估计错误。MUSIC谱图随快拍数变化情况在上面的实验中，快拍数为5000,分辨效果较好，原因是当快拍足够时，能 较好地得到噪声的统计特性，从而精确地分解出信号子空间和噪声子空间， 从而 较好的分辨信源方向，当快拍数下降时，分辨效果将会下降。-5-10-15-20-25-30-100-80-60-40-20020406080100快拍数为300快拍数为100可见当快拍数下降的时候，峰值将会没有那么尖锐，对信号的分辨力下降， 甚至无法两个角度差较小的信号，原因是小样本时间平均与统计平均差距较大， 噪声子空间分解不准确。5、实验结论由上面的实验可以看出，music算法可以实现对于DO

9、A的求解，信噪比越大,快拍数越高越有利于分辨力的提高。 此外不同的music算法求解时，计算时间和 精度都有差别，但是由于来波来自3个方向，每个方向的求解实际上是存在相互 牵制的关系，比如来自0的信号求解很准可能就会造成-20。的信号求解不准。 而信号的窄带信也是很关键的，当不满足窄带条件时，求解的角度会误差很大。6、实验程序 clear all ;clc;close M=16;K=3; %D? e yNN=M-K; %?! u e y m=0:1:M-1;wavelength=1; %f=300MHzf=3*10A8/wavelength;fs=10*f;Ts=1/fs;df=f/30; %

10、df=10MHz d=wavelength/2;%?5 a D?ak=2*pi/wavelength;%2，e yN=500; %?i ?e y%?-?-?-?- e e o ?D?o?1?a?Q X ?-?-?-?-%U?1|? e e /dBU?1|? e e /dB %?登 u 1|? eSNR=-20; %D? e dB SCR2=0; %D? 2o ?D?7 SCR3=0; %D? 3o ?D? noise=10A(-SNR/10);pj2=10A(SCR2/10); pj3=10A(SCR3/10); as2=sqrt(pj2);as3=sqrt(pj3); sigma=sqrt(

11、noise); fai_j1=2/180*pi;fai_j2=30/180*pi;fai_j3=31/180*pi;w1=f*2*pi;w2=(f+df)*2*pi;w3=(f-df)*2*pi;kj1=k*sin(fai_j1);%D? 21|? e%D? 31|? e%D? 2 - u ?e%D? 3 - u ?e%?! u , u ? e%D? ? 1?e%D? 2?e%D? 3?ekj2=k*w2/w1*sin(fai_j2);kj3=k*w3/w1*sin(fai_j3);akj1=exp(-1j*(kj1*m*d).;akj2=exp(-1j*(kj2*m*d).;akj3=exp

12、(-1j*(kj3*m*d).;R=zeros(M,M);Ri=zeros(M,M);for i=1:Nri=normrnd(0,sigma,M,1);Ri=(akj1*exp(1j*w1*i*Ts)+as2*akj2*exp(1j*w2*i*Ts)+as3*akj3*exp(1j*w3*i*Ts)+ri)*(akj1*exp(1j*w1*i*Ts)+as2*akj2*exp(1j*w2*i*Ts)+as3*akj3*exp(1j*w3*i*Ts)+ri);R=R+Ri;endR=R/N;V,D=eigs(R,NN, sm);fai=(-90:0.005:90)/180*pi;ps=zeros

13、(1,length(fai);disp( ?- Cmusic? . ?a 2 ?e /? e);for i=1:length(fai)%?6 a Da + D?e ？ 6ki=k*w3/w1*sin(fai(i);aki=exp(-1j*(ki*m*d).;ps(i)=real(1/(aki*(V*V)*aki);endfor i=2:length(fai)-1%? ?o ? 6?小仙？ ? a 2 ?oif ps(i)ps(i+1)&ps(i)ps(i-1)&ps(i)max(ps)/20 disp(fai(i)/pi*180);泡2. ?6endendfigure;plot(fai/pi*

14、180,10*log10(ps/max(ps);%? u ? %syms z ;Pz1=sym(zeros(1,M);Pz2=sym(zeros(1,M).;for i=1:MPz1(i)=zA(1-i);Pz2(i)=zA(i-1); endtsym=solve(Pz1*(V*V)*Pz2);tdouble=double(tsym);r=abs(tdouble);rnum=length(r);%numsolve=zeros(1,K);% for i=1:K%? 0 ?o o?e e ? ?u%temp1=20;% temp2=1;% for j=1:length(r)%if abs(r(j)

15、-1)temp1%temp1=abs(r(j)-1);%temp2=j;%end% % %endendnumsolve(i)=tdouble(temp2); r(temp2)=10;tsolve=zeros(1,rnum);j=0;for i=1:rnumif (abs(r(i)-1)0.1) j=j+1;tsolve(j)=tdouble(i);endendnumsolve=tsolve(1:j);signal_position=asin(-wavelength/(2*pi*d)*angle(numsolve) )/pi*180;real_position=signal_position(1);for