MUSIC算法频率估计

采用MUSIC方法的白噪声频率检测仿真本试验提供了一种使用MUSIC方法的白噪声中一个正弦信号和M 个正弦信号的特征分解频率估计的仿真试验，并讨论了虚假峰的成因并给出了实验证明。

问题描述假定仿真的观测数据分别由 （1）单个正弦信号检测的情况()43()4()j n x n eu n ππ+=+（2）多个正弦信号检测的情况5()()()433645()423()j n j n j n x n eeeu n ππππππ+++=+++产生，其中是一高斯白噪声，其均值为0，方差为1。

用MUSIC 方法估计观测数据中正弦波的频率，并给出白噪声方差()u n 2u σ 与复正弦波的振幅A 的估计值。

多重信号分类的MUSIC 方法实际应用中常常需要对空间中存在的多个信号源进行分解，以便跟踪或检测我们感兴趣的空间信号，抑制那些被认为是干扰的空间信号。

对天线阵列接收的空间信号所进行的分析与处理称为阵列信号处理。

而空间谱估计技术是在波束形成技术、零点技术和时域谱估计技术的基础上发展起来的一种技术。

与频谱表示信号在各个频率上的能量分布相对应，空间谱则可解释为信号在空间各个方向上的能量分布，空间谱估计技术的目标是研究提高在处理带宽内空间信号角度的估计精度、角度分辨率和提高运算速度的各种算法。

经过多年的发展，已经产生了大量性能优异的测向算法可资利用，典型的有MUSIC.ESPRIT、子空间拟合、多维MUSIC 等。

MUSIC 算法是基于特征结构分析的空间谱估计方法，是空间谱估计技术的典型代表。

其测向原理是根据矩阵特征分解的理论，对阵列输出协方差矩阵进行特征分解，将信号空间分解为噪声子空间G 和信号子空间S，利用噪声子空间G 与阵列的方向矩阵A 的列矢量正交的性质，构造空间谱函数P(w)并进行谱峰搜索，从而估计出波达方向信息。

设空间有p 个互不相关的信号以方位角12,,p θθ""θ入射到具有m 个接收阵元的接收阵元阵列中，入射信号的数目p 小于阵列的阵元数m。

则此阵列系统的信号模型为：1()()()()()()()pi i i x n a w s n u n A w s n u ==+=∑n +其中12121(1)(1)(1)()[(),,()]111p p p jwjw jw j m w j m w j m w A w a w a w ee e e e e −−−−−−−−−=⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦"""####" 1()[(),,()]p s n s n s n ="对由上式描述的阵列信号观测模型做以下假设： 假设1：对于不同的值，向量相互线性独立；i w ()i a w 假设2：加性噪声向量的每个原色都是零均值的复白噪声，它们不相关，并且具有相同的方差()u n 2u σ；假设3：矩阵非奇异，即{()()}H P E s n s n =()rank P p =。

上述三个假设都只是一般的假设，在实际中容易得到满足。

可见，22{()()}(){()()}() H xx H H u H u R E x n x n A w E s n s n A w I APA Iσσ==+=+ 于是xx R 是个对称矩阵，令其特征值分解为H xx R U U =Σ式中221(,,m diag )σσΣ="。

由于A 满列秩，故()()H rank APA rank P p ==，这里假定p m <。

于是有221(,,,0,,0)H H p U APA U diag αα=""其中21,,p2αα"是无加性噪声时的观测信号()Ax n 的自相关矩阵H APA 的特征值。

同时用左乘和U 右乘上式，得H U 2221 (,,,0,,0)H H H H xx u p u U R U U APA U U Udiag Iσ2αασ=+=+""这表明，自相关矩阵xx R 的特征值为2222, i=1,, i=1,,i u i iu pp mασλσσ⎧+⎪==⎨+⎪⎩"" 根据信号特征值和噪声特征值，将特征矩阵U 的列向量分成两部分，即[]U S G =#其中1111[,,][,,][,,][,,]p p m p p m S s s u u G g g u u −+====""""其中分别由信号特征向量和噪声特征向量组成。

考察2[][]H xx u H S O R G S G S G G I G σ⎡⎤⎡⎤=Σ=Σ=⎢⎥⎢⎥⎣⎦⎣⎦##又由2H xx u R APA I σ=+有2H xx u R G APA G G σ=+，利用上式的结果，得到H APA G O =进而有H H G APA G O =上式成立的充分必要条件是H A G O =将1[(),,()]p A a w a w ="要代入上式得1()0,,,H T p w G w w w α=="显然，当时，1,,p w w w ≠"()0H T w G α≠。

将上式改成标量形式，可以定义一种类似于功率谱的函数：1()()()H H P w w GG w αα=上式取峰值的p 个w 值给出p 个信号的波达方向12,,p θθθ""。

这样定义的函数描述了空间参数（即波达方向）的分布，因此成为空间谱。

由于它能够对多个信号进行识别，故以此得名MUSIC 方法（multiple signal classification）。

白噪声中单个正弦信号的频率检测与估计1． 实验结果（1）短数据时的估计结果，N_x=64.w1_estimate = 0.2500（π） （估计频率）S = 1.4058 （噪声功率）A = 3.9835 （信号振幅估计）（2）长数据时的估计结果，N_x=1024。

w1_estimate = 0.2500（π） （估计频率）S = 0.9646 （噪声功率）A = 4.0083 （信号振幅估计）2．对结果的讨论（1）由给出的试验结果可知，在π/4处有明显的峰值，其它地方都是平坦的，频率估计效果非常好。

（2）在信号长度很短（64）的情况下，σ2 和A 的估计误差比较大。

造成误差大的原因是：由于信号短，加窗效应较明显，自相关矩阵的估计误差比较大。

增加信号长度可以减小自相关矩阵的估计误差，可以得到比较理想的结果，这可以从（N_x=1024）的结果中看出。

白噪声中多个正弦信号的频率检测与估计1．实验结果（1）短数据时的估计结果，N_x=64.S = 2.1681 （噪声功率）A = 3.8872 2.1874 2.9104 （信号振幅估计）（2）短数据时的估计结果，N_x=1024.S = 1.0385 （噪声功率）A = 4.0098 2.0218 2.9917 （信号振幅估计） 2．对结果的讨论（1）由给出的试验结果可知，在π/4，π/3，5*π/4 处分别都有明显的峰，与实际值十分吻合。

为了节省篇幅没有计算估计值，只从图上观察是否相符，有兴趣的可以分段找最大值的方法将估计值算出来。

（2）由信号长度导致的噪声功率信号振幅估计的差别在上述结果中也有明显的体现，产生的原因同上述分析。

（3）疑问：从图中发现三个峰高度有明显的差异，开始以为与Ai 或ψi 有关，经过实验并未发现有何联系，未找到影响峰值的原因。

虚假峰成因与数量的讨论由于在MUSIC方法中，前面的几个噪声特征值对应的特征向量都是相同的，于是尝试用单个特征矢量进行频率估计。

在用单个特征矢量估计ω1 的时候会出现虚假峰，即不等于ω1 的位置上多出几个峰。

研究特征向量的正交性可知，特征矢量是两两正交的。

下面的程序可以证实：对特征分解得到的特征矢量组成的矩阵V 进行处理，判断列矢量之间，行矢量之间是否正交：% Test_V.m查看特征矢量是否互相垂直。

结果为：z=zeros(6); z = 1 0 0 0 0 0for i=1:6 0 1 0 0 0 0 for j=i:6 0 0 1 0 0 0z(i,j)=V(:,i)'*V(:,j); 0 0 0 1 0 0z(j,i)=V(j,:)*V(i,:)'; 0 0 0 0 1 0 end 0 0 0 0 0 1end。

这说明特征矢量矩阵的行列矢量（即特征矢量）之间，行矢量之间都是两两正交的。

特征矢量两两正交，也就是说eHVi 应该有N-1 个零点，因此应该有N-1 个峰值，这是虚假峰产生的真正原因！为此特作了多次实验加以验证，模型仍用前面1 个复正弦信号，ω1=π/4，N=6，结果如下图所示：通过上述试验加深了对MUSIC方法的认识。

结果显示说明mymusics.m 白噪声中单个正弦信号的频率检测与估计clear;close all;%Frequency Estimation by Eigendecomposition of Autocorrelation Matrix N_x=1024; % Length of SignalN=10; % Size of Rx MatrixA=[4 2 3];w=[pi/4 pi/3 5*pi/4]';phase=[pi/3*ones(1,N_x);pi/6*ones(1,N_x);pi/5*ones(1,N_x)];M=3; % Number of Signalsx=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([0 2*pi])set(gca,'XTick',0:pi/4:2*pi)set(gca,'XTickLabel',{'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;close all;%Frequency Estimation by Eigendecomposition of Autocorrelation Matrix N_x=1024; % Length of SignalN=6; % Size of Rx Matrixw1=pi/4;x=4*exp(j*(w1*[0:N_x-1]+pi/3))+randn(1,N_x); %Generate the signal: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([0 2*pi])set(gca,'XTick',0:pi/4:2*pi)set(gca,'XTickLabel',{'0','pi/4','pi/2','3pi/4','pi','5pi/4','3pi/2' ,'7pi/4','2pi'})[max I]=max(Pw); %find the index of max Pww1_estimate=(I-1)/Nw %compute the estimate of w1 by the indexS=mean(D(1:N-1))A=sqrt(D(6)/N)致谢本报告参考了网上张建军的实践报告，并与其进行了简单的讨论，感谢他的工作和讨论给出的帮助！。