MUSIC算法频率估计

1、采用MUSIC方法的白噪声频率检测仿真本试验提供了一种使用MUSIC方法的白噪声中一个正弦信号和M个正弦信号的特征分解频率估计的仿真试验，并讨论了虚假峰的成因并给出了实验证明。问题描述假定仿真的观测数据分别由单个正弦信号检测的情况vnn、x(n)=4ej(4n+3)+u(n)多个正弦信号检测的情况nnnn5nn、x(n)=4ej(4n+3)+2ej(3n+6)+3ej(4n+5)+u(n)产生，其中u(n)是一高斯白噪声，其均值为0,方差为1。用MUSIC方法估计观测数据中正弦波的频率，并给出白噪声方差Q2与复正弦波的振幅A的估计值。u多重信号分类的MUSIC方法实际应用中常常需要对空间中存在

2、的多个信号源进行分解，以便跟踪或检测我们感兴趣的空间信号，抑制那些被认为是干扰的空间信号。对天线阵列接收的空间信号所进行的分析与处理称为阵列信号处理。而空间谱估计技术是在波束形成技术、零点技术和时域谱估计技术的基础上发展起来的一种技术。与频谱表示信号在各个频率上的能量分布相对应，空间谱则可解释为信号在空间各个方向上的能量分布，空间谱估计技术的目标是研究提高在处理带宽内空间信号角度的估计精度、角度分辨率和提高运算速度的各种算法。经过多年的发展，已经产生了大量性能优异的测向算法可资利用，典型的有MUSIC.ESPRIT、子空间拟合、多维MUSIC等。MUSIC算法是基于特征结构分析的空间谱估计方法

3、，是空间谱估计技术的典型代表。其测向原理是根据矩阵特征分解的理论，对阵列输出协方差矩阵进行特征分解，将信号空间分解为噪声子空间G和信号子空间S,利用噪声子空间G与阵列的方向矩阵A的列矢量正交的性质，构造空间谱函数P(w)并进行谱峰搜索，从而估计出波达方向信息。设空间有p个互不相关的信号以方位角9,9,9入射到具有m个接收阵元12p的接收阵元阵列中，入射信号的数目p小于阵列的阵元数m。则此阵列系统的信号模型为：x(n)=区a(w)s(n)+u(n)=A(w)s(n)+u(n)iii=1其中A(w)=a(w),a(w)1p_111_e-jWe-丿叫e-jwpe-j(m-i)We-j(m-i)w?e

4、j(m-1)wps(n)=s(n),s(n)1p对由上式描述的阵列信号观测模型做以下假设：假设1：对于不同的w值，向量a(w)相互线性独立；ii假设2：加性噪声向量u(n)的每个原色都是零均值的复白噪声，它们不相关，并且具有相同的方差a2；u假设3：矩阵非奇异P=Es(n)sh(n)，即rank(P)=p。上述三个假设都只是一般的假设，在实际中容易得到满足。可见，R=Ex(n)xH(n)xx=A(w)Es(n)sH(n)AH(w)+a2Iu=APAH+a2Iu于是R是个对称矩阵，令其特征值分解为xxR=UZUhxx式中Z=diag(a2,.,a2)。1m由于A满列秩，故rank(APAh)=r

5、ank(P)=p，这里假定pm。于是有UhAPAhU=diag(a2,a2,0,0)1p其中a2,a2是无加性噪声时的观测信号Ax(n)的自相关矩阵APAh的特征值。1p同时用Uh左乘和U右乘上式，得UhRU=UhAPAhU+a2UhUxxu=diag(a2,a2,0，,0)+a211pu这表明，自相关矩阵R的特征值为xxra2+a2,i=1,p九=O2=iuio2i=p+1,mu根据信号特征值和噪声特征值，将特征矩阵U的列向量分成两部分，即U=S#G其中S=s,s=u,u1p1pG=g,g=u,u1m-pp+1m其中分别由信号特征向量和噪声特征向量组成。考察RG=S#GZxxSHGH=O2G

6、u又由R=APAh+o21有RG=APAhG+o2G,利用上式的结果，得到xxuxxuAPAHG=O进而有GHAPAHG=O上式成立的充分必要条件是AHG=O将A=a(w),a(w)要代入上式得1paH(w)G=0T,w=w,w1p显然，当w丰w，,w时，ah(w)G丰0t。1p将上式改成标量形式，可以定义一种类似于功率谱的函数1ah(w)GGha(w)上式取峰值的p个w值给出p个信号的波达方向9,9,9。12p这样定义的函数描述了空间参数(即波达方向)的分布，因此成为空间谱。由于它能够对多个信号进行识别，故以此得名MUSIC方法(multiplesignalclassification)。白

7、噪声中单个正弦信号的频率检测与估计1实验结果(1)短数据时的估计结果，N_x=64.S=1.4058A=3.9835(2)长数据时的估计结果，N_x=1024。噪声功率)信号振幅估计)wl_estimate=0.2500(n)(估计频率)0.9646(噪声功率)00pi/42对结果的讨论pi/23pi/4(信号振幅估计)5pi/43pi/27pr/42pi(1)由给出的试验结果可知，在n/4处有明显的峰值，其它地方都是平坦的，频率估计效果非常好。(2)在信号长度很短(64)的情况下，a2和A的估计误差比较大。造成误差大的原因是：由于信号短，加窗效应较明显，自相关矩阵的估计误差比较大。增加信号长

8、度可以减小自相关矩阵的估计误差，可以得到比较理想的结果，这可以从(N_x=1024)的结果中看出。5pi/100白噪声中多个正弦信号的频率检测与估计1实验结果（1）短数据时的估计结果，N_x=64.700600500400（噪声功率）300-S=2.1681(2)A=3.88722.18742.910短数据时的估计结果，N_x=1（24.（信号振幅估计）0pi/4pi/23pi/4Pi5p80007000600050004000S=1.03853000-A=4.0098(噪声功率)2.02182.9917(信号振幅估计)2对结果的讨论(1)由给出的试验结果可知，在n/4,n/3,5*n/4处分

9、别都有明显的峰,与实际值十分吻合。为了予省篇幅没有趣的可以分段找最大值的方法将估计值与实际值十分吻合计算估计值，只从图上观察是否相符，有兴5pi/1005pi/100由信号长度导致的噪声功率信号振幅估计的差别在上述结果中也有明显的体现，产生的原因同11L0pi/4pi/23pi/4pi疑问：从图中发现三个峰高度有明显的差异，开始以为与Ai或有关，经过实验并未发现有何联系，未找到影响峰值的原因。100虚假峰成因与数量的讨论由于在MUSIC方法中，前面的几个噪声特征值对应的特征向量都是相同的，于是尝试用单个特征矢量进行频率估计。在用单个特征矢量估计31的时候会出现虚假峰，即不等于31的位置上多出几

10、个峰。研究特征向量的正交性可知，特征矢量是两两正交的。下面的程序可以证实：对特征分解得到的特征矢量组成的矩阵V进行处理，判断列矢量之间，行矢量之间是否正交：%Test_V.m查看特征矢量是否互相垂直。z=zeros(6);结果为：z=100000fori=1:6010000forj=i:6z(i,j)=V(:,i)*V(:,j);z(j,i)=V(j,:)*V(i,:);001000000100000010end000001end。这说明特征矢量矩阵的行列矢量(即特征矢量)之间，行矢量之间都是两两正交的。特征矢量两两正交，也就是说eHVi应该有N-1个零点，因此应该有N-1个峰值，这是虚假峰产

11、生的真正原因！为此特作了多次实验加以验证，模型仍用前面1个复正弦信号，31=n/4，N=6，结果如下图所示：mymusics.m白噪声中单个正弦信号的频率检测与估计clear;closeall;%FrequencyEstimationbyEigendecompositionofAutocorrelationMatrixN_x=1024;%LengthofSignalN=10;%SizeofRxMatrixA=423;w=pi/4pi/35*pi/4;phase=pi/3*ones(1,N_x);pi/6*ones(1,N_x);pi/5*ones(1,N_x);M=3;%NumberofSig

12、nalsx=randn(1,N_x)+A*exp(j*(w*0:N_x-1+phase);Cx=xcorr(x,biased);Rxx=Cx(N_x:N_x+N-1);Rx=toeplitz(Rxx);V,D=eig(Rx);%Eigendecomposition特征分解D=sum(D);Nw=128;ww=0:256/128*pi;e=exp(-j*ww*0:N-1);ev=e*V(:,1:N-M);Pw=1./real(diag(ev*ev);figureplot(ww,Pw);xlim(02*pi)set(gca,XTick,0:pi/4:2*pi)set(gca,XTickLabel,

13、0,pi/4,pi/2,3pi/4,pi,5pi/4,3pi/2,7pi/4,2pi)S=mean(D(1:N-M)E=exp(j*0:N-1.*w);P=real(E(Rx-(eye(N).*S)/E);A=sqrt(sum(P)mymusicm.m白噪声中多个正弦信号的频率检测与估计clear;closeall;%FrequencyEstimationbyEigendecompositionofAutocorrelationMatrixN_x=1024;%LengthofSignalN=6;%SizeofRxMatrixw1=pi/4;x=4*exp(j*(w1*0:N_x-1+pi/3)+randn(1,N_x);%Generatethesignal:Cx=xcorr(x,biased);Rxx=Cx(N_x:N_x+N-1);Rx=toeplitz(Rxx);V,D=eig(Rx);%Eigendecomposition特征分解D=sum(D);Nw=128;w=0:256/128*pi;e=exp(-j*w*0:N-1);ev=e*V(:,1:N-1);Pw=1./real(diag(ev*ev);plot(w,Pw);xlim(02*pi)set(gca,XTick,0:pi/4:2*pi)se