较为详细的MUSIC算法原理及MATLAB实现

阵列信号处理 专栏收录该内容
2 篇文章 1 订阅

DOA估计算法

DOA(Direction Of Arrival)波达方向定位技术主要有ARMA谱分析、最大似然法、熵谱分析法和特征分解法,特征分解法主要有MUSIC算法、ESPRIT算法WSF算法等。
MUSIC (Multiple Signal Classification)算法,即多信号分类算法,由Schmidt等人于1979年提出。MUSIC算法是一种基于子空间分解的算法,它利用信号子空间和噪声子空间的正交性,构建空间谱函数,通过谱峰搜索,估计信号的参数。对于声源定位来说,需要估计信号的DOA。MUSIC算法对DOA的估计有很高的分辨率,且对麦克风阵列的形状没有特殊要求,因此应用十分广泛。
在这里插入图片描述

运用矩阵的定义,可得到更为简洁的表达式:

X = A S + N X=AS+N X=AS+N

式中

X = [ x 1 ( t ) , x 2 ( t ) , . . . x M ( t ) ] T X=[x_1(t),x_2(t),...x_M(t)]^T X=[x1(t),x2(t),...xM(t)]T
S = [ S 1 ( t ) , S 2 ( t ) , . . . S D ( t ) ] T S=[S_1(t),S_2(t),...S_D(t)]^T S=[S1(t),S2(t),...SD(t)]T
A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T A=[a(\theta_1),a(\theta_2),...a(\theta_D)]^T A=[a(θ1),a(θ2),...a(θD)]T
N = [ n 1 ( t ) , n 2 ( t ) , . . . n M ( t ) ] T N=[n_1(t),n_2(t),...n_M(t)]^T N=[n1(t),n2(t),...nM(t)]T

X X X为阵元的输出, A A A为方向响应向量, S S S是入射信号, N N N表示阵列噪声。

其中 φ k = 2 π d λ s i n θ k \varphi_k=\frac{2\pi d}{\lambda}sin\theta_k φk=λ2πdsinθk

A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T = [ 1 1 ⋯ 1 e − j φ 1 e − j φ 2 ⋯ e − j φ D ⋮ ⋮ ⋱ ⋮ e − j ( M − 1 ) φ 1 e − j ( M − 1 ) φ 2 ⋯ e − j ( M − 1 ) φ D ] A=[a(\theta_1),a(\theta_2),...a(\theta_D)]^T=\left[ \begin{matrix} 1 & 1 & \cdots & 1 \\ {e^{-j\varphi_1}} & {e^{-j\varphi_2}} & \cdots & {e^{-j\varphi_D}} \\ \vdots & \vdots & \ddots & \vdots \\ {e^{-j(M-1)\varphi_1}} & {e^{-j(M-1)\varphi_2}} & \cdots & {e^{-j(M-1)\varphi_D}} \\ \end{matrix} \right] A=[a(θ1),a(θ2),...a(θD)]T=1ejφ1ej(M1)φ11ejφ2ej(M1)φ21ejφDej(M1)φD
x m ( t ) x_m(t) xm(t)进行N点采样,要处理的问题就变成了通过输出信号 x m ( t ) x_m(t) xm(t)的采样 { x m ( i ) = 1 , 2 , . . . , M } \{ x_m (i)=1,2,...,M\} {xm(i)=1,2,...,M}估计信号源的波达方向角 θ 1 , θ 2 . . . θ D \theta_1,\theta_2...\theta_D θ1,θ2...θD,由此可以很自然的将阵列信号看作是噪声干扰的若干空间谐波的叠加,从而将波达方向估计问题与谱估计联系起来。
对阵列输出X做相关处理,得到其协方差矩阵

R x = E [ X X H ] R_x=E[XX^H] Rx=E[XXH]

其中 H H H表示矩阵的共轭转置。
根据已假设信号与噪声互不相关、噪声为零均值白噪声,因此可得到:

R x = E [ ( A S + N ) ( A S + N ) H ] = A E [ S S H ] A H + E [ N N H ] = A R S A H + R N R_x=E[(AS+N)(AS+N)^H] =AE[SS^H]A^H+E[NN^H]=AR_SA^H+R_N Rx=E[(AS+N)(AS+N)H]=AE[SSH]AH+E[NNH]=ARSAH+RN
其中 R s = E [ S S H ] R_s=E[SS^H] Rs=E[SSH]称为信号相关矩阵

R N = σ 2 I R_N=\sigma^2I RN=σ2I是噪声相关阵

σ 2 \sigma^2 σ2是噪声功率

I I I M × M M\times M M×M阶的单位矩阵

在实际应用中通常无法直接得到 R x R_x Rx,能使用的只有样本的协方差矩阵:
R x ^ = 1 N ∑ i = 1 N X ( i ) X H ( i ) \hat{R_x}=\frac{1}{N} \sum_{i=1}^{N}X(i)X^H (i) Rx^=N1i=1NX(i)XH(i) R x ^ \hat{R_x} Rx^ R x R_x Rx的最大似然估计。

当采样数 N → ∞ N\to\infty N,他们是一致的,但实际情况将由于样本数有限而造成误差。根据矩阵特征分解的理论,可对阵列协方差矩阵进行特征分解,首先考虑理想情况,即无噪声的情况: R x = A R s A H R_x=AR_sA^H Rx=ARsAH,对均匀线阵,矩阵A由
A = [ a ( θ 1 ) , a ( θ 2 ) , . . . a ( θ D ) ] T = [ 1 1 ⋯ 1 e − j φ 1 e − j φ 2 ⋯ e − j φ D ⋮ ⋮ ⋱ ⋮ e − j ( M − 1 ) φ 1 e − j ( M − 1 ) φ 2 ⋯ e − j ( M − 1 ) φ D ] A=[a(\theta_1),a(\theta_2),...a(\theta_D)]^T=\left[ \begin{matrix} 1 & 1 & \cdots & 1 \\ {e^{-j\varphi_1}} & {e^{-j\varphi_2}} & \cdots & {e^{-j\varphi_D}} \\ \vdots & \vdots & \ddots & \vdots \\ {e^{-j(M-1)\varphi_1}} & {e^{-j(M-1)\varphi_2}} & \cdots & {e^{-j(M-1)\varphi_D}} \\ \end{matrix} \right] A=[a(θ1),a(θ2),...a(θD)]T=1ejφ1ej(M1)φ11ejφ2ej(M1)φ21ejφDej(M1)φD
所定义的范德蒙德矩阵,只要满足 θ i ≠ θ j , i ≠ j \theta_i\neq \theta_j,i\neq j θi=θj,i=j,则他的各列相互独立。

R s R_s Rs为非奇异矩阵 R a n k ( R s ) = D Rank(R_s)=D Rank(Rs)=D,各信号源两两不相干,且 M > D M>D M>D,则 r a n d ( A R s A H ) = D rand(AR_sA^H)=D rand(ARsAH)=D

由于 R x = E [ X X H ] R_x=E[XX^H] Rx=E[XXH],有:

R s H = R x R_s^H=R_x RsH=Rx

R s R_s Rs为Hermite矩阵,它的特性是都是实数,又由于 R s R_s Rs为正定的,因此 A R s A … … H AR_sA……H ARsAH为半正定的,它有D个正特征值和 M − D M-D MD个零特征值。

再考虑有噪声存在的情况

R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2I Rx=ARsAH+σ2I

由于 σ 2 > 0 \sigma^2>0 σ2>0 R x R_x Rx为满秩阵,所以 R x R_x Rx有M个正实特征值 λ 1 , λ 2 . . . λ M \lambda_1,\lambda_2...\lambda_M λ1,λ2...λM

分别对应于M个特征向量 v 1 , v 2 . . . v M v_1,v_2...v_M v1,v2...vM。又由于 R x R_x Rx为Hermite矩阵,所以各特征向量是正交的,即: v i H v j = 0 , i ≠ j v_i^Hv_j=0,i\neq j viHvj=0,i=j与信号有关的特征值只有D个,分别等于矩阵 A R s A H AR_sA^H ARsAH的各特征值与 σ 2 \sigma^2 σ2之和,其余 M − D M-D MD个特征值为 σ 2 \sigma^2 σ2,即 σ 2 \sigma^2 σ2 R R R的最小特征值,它是 M − D M-D MD维的,对应的特征向量 v i , i = 1 , 2 , . . . , M v_i,i=1,2,...,M vi,i=1,2,...,M中,也有D个是与信号有关的,另外 M − D M-D MD个是与噪声有关的,可利用特征分解的性质求出信号源的波达方向 θ k \theta_k θk

MUSIC算法的原理及实现

通过对协方差矩阵的特征值分解,可得到如下结论:

将矩阵 R x R_x Rx的特征值进行从小到大的排序,即 λ 1 ≥ λ 2 ≥ . . . ≥ λ M > 0 \lambda_1 \geq \lambda_2\geq...\geq\lambda_M>0 λ1λ2...λM>0,其中D个较大的特征值对应于信号, M − D M-D MD个较小的特征值对应于噪声。

矩阵 R x R_x Rx的属于这些特征值的特征向量也分别对应于各个信号和噪声,因此可把 R x R_x Rx的特征值(特征向量)划分为信号特征(特征向量)与噪声特征(特征向量)。

λ i \lambda_i λi R x R_x Rx的第 i i i个特征值, v i v_i vi是与 λ i \lambda_i λi个相对应的特征向量,有:

R x v i = λ i v i R_xv_i=\lambda_iv_i Rxvi=λivi

再设 λ i = σ 2 \lambda_i=\sigma^2 λi=σ2 R x R_x Rx的最小特征值 R x v i = σ 2 v i i = D + 1 , D + 2... M R_xv_i=\sigma^2v_i i=D+1,D+2...M Rxvi=σ2vii=D+1,D+2...M

R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2I Rx=ARsAH+σ2I代入可得 σ 2 v i = ( A R s A H + σ 2 I ) v i \sigma^2v_i=(AR_sA^H+\sigma^2I)v_i σ2vi=(ARsAH+σ2I)vi

将其右边展开与左边比较得:

A R s A H v i = 0 AR_sA^Hv_i=0 ARsAHvi=0

A H A A^HA AHA D ∗ D D*D DD维的满秩矩阵, ( A H A ) − 1 (A^HA)^{-1} (AHA)1存在;

R s − 1 R_s^{-1} Rs1同样存在,则上式两边同乘以 R s − 1 ( A H A ) − 1 A H R_s^{-1}(A^HA)^{-1}A^H Rs1(AHA)1AH

有:

R s − 1 ( A H A ) − 1 A H A R s A H v i = 0 R_s^{-1}(A^HA)^{-1}A^HAR_sA^Hv_i=0 Rs1(AHA)1AHARsAHvi=0

于是有

A H v i = 0 , i = D + 1 , D + 2 , . . . , M A^Hv_i=0,i=D+1,D+2,...,M AHvi=0i=D+1,D+2,...,M

上式表明:噪声特征值所对应的特征向量(称为噪声特征向量) v i v_i vi,与矩阵 A A A的列向量正交,而 A A A的各列是与信号源的方向相对应的,这就是利用噪声特征向量求解信号源方向的出发点。

用各噪声特征向量为例,构造一个噪声矩阵 E n E_n En:

E n = [ v D + 1 , v D + 2 , . . . v M ] E_n=[v_{D+1},v_{D+2},...v_{M}] En=[vD+1,vD+2,...vM]

定义空间谱 P m u ( θ ) P_{mu}(\theta) Pmu(θ):

P m u ( θ ) = 1 a H ( θ ) E n E n H a ( θ ) = 1 ∥ E n H a ( θ ) ∥ 2 P_{mu}(\theta)=\frac{1}{{a^H}(\theta)}E_nE_n^Ha(\theta)=\frac{1}{\Vert E_n^Ha(\theta)\Vert^2} Pmu(θ)=aH(θ)1EnEnHa(θ)=EnHa(θ)21

该式中分母是信号向量和噪声矩阵的内积,当 a ( θ ) a(\theta) a(θ) E n E_n En的各列正交时,该分母为零,但由于噪声的存在,它实际上为一最小值,因此 P m u ( θ ) P_{mu}(\theta) Pmu(θ)有一尖峰值,由该式,使 θ \theta θ变化,通过寻找波峰来估计到达角。

MUSIC算法实现的步骤

1.根据N个接收信号矢量得到下面协方差矩阵的估计值:

R x = 1 N ∑ i = 1 N X ( i ) X H ( i ) R_x=\frac{1}{N} \sum_{i=1}^NX(i)X^H(i) Rx=N1i=1NX(i)XH(i)

对上面得到的协方差矩阵进行特征分解

R x = A R s A H + σ 2 I R_x=AR_sA^H+\sigma^2I Rx=ARsAH+σ2I

2.按特征值的大小排序将与信号个数 D D D相等的特征值和对应的特征向量看做信号部分空间,将剩下的 M − D M-D MD个特征值和特征向量看做噪声部分空间,得到噪声矩阵 E n E_n En:

A H v i = 0 , i = D + 1 , D + 2 , . . . , M A^Hv_i=0,i=D+1,D+2,...,M AHvi=0i=D+1,D+2,...,M

E n = [ v D + 1 , v D + 2 , . . . v M ] E_n=[v_{D+1},v_{D+2},...v_{M}] En=[vD+1,vD+2,...vM]

3.使 θ \theta θ变化,按照式

P m u ( θ ) = 1 a H ( θ ) E n E n H a ( θ ) P_{mu}(\theta)=\frac{1}{{a^H}(\theta)E_nE_n^Ha(\theta)} Pmu(θ)=aH(θ)EnEnHa(θ)1

来计算谱函数,通过寻求峰值来得到波达方向的估计值。

clear; close all;
%%%%%%%% MUSIC for Uniform Linear Array%%%%%%%%
derad = pi/180;      %角度->弧度
N = 8;               % 阵元个数        
M = 3;               % 信源数目
theta = [-30 0 60];  % 待估计角度
snr = 10;            % 信噪比
K = 512;             % 快拍数
 
dd = 0.5;            % 阵元间距 
d=0:dd:(N-1)*dd;
A=exp(-1i*2*pi*d.'*sin(theta*derad));  %方向矢量

%%%%构建信号模型%%%%%
S=randn(M,K);             %信源信号,入射信号
X=A*S;                    %构造接收信号
X1=awgn(X,snr,'measured'); %将白色高斯噪声添加到信号中
% 计算协方差矩阵
Rxx=X1*X1'/K;
% 特征值分解
[EV,D]=eig(Rxx);                   %特征值分解
EVA=diag(D)';                      %将特征值矩阵对角线提取并转为一行
[EVA,I]=sort(EVA);                 %将特征值排序 从小到大
EV=fliplr(EV(:,I));                % 对应特征矢量排序
                 
 
% 遍历每个角度,计算空间谱
for iang = 1:361
    angle(iang)=(iang-181)/2;
    phim=derad*angle(iang);
    a=exp(-1i*2*pi*d*sin(phim)).'; 
    En=EV(:,M+1:N);                   % 取矩阵的第M+1到N列组成噪声子空间
    Pmusic(iang)=1/(a'*En*En'*a);
end
Pmusic=abs(Pmusic);
Pmmax=max(Pmusic)
Pmusic=10*log10(Pmusic/Pmmax);            % 归一化处理
h=plot(angle,Pmusic);
set(h,'Linewidth',2);
xlabel('入射角/(degree)');
ylabel('空间谱/(dB)');
set(gca, 'XTick',[-90:30:90]);
grid on;

实现结果

在这里插入图片描述

©️2021 CSDN 皮肤主题: 技术黑板 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值