DoA 估计:多重信号分类 MUSIC 算法(附 MATLAB 代码)

这篇具有很好参考价值的文章主要介绍了DoA 估计:多重信号分类 MUSIC 算法(附 MATLAB 代码)。希望对大家有所帮助。如果存在错误或未考虑完全的地方,请大家不吝赐教,您也可以点击"举报违法"按钮提交疑问。

本文首次在公众号【零妖阁】上发表,为了方便阅读和分享,我们将在其他平台进行自动同步。由于不同平台的排版格式可能存在差异,为了避免影响阅读体验,建议如有排版问题,可前往公众号查看原文。感谢您的阅读和支持!

DoA 估计是指根据天线阵列的接收信号估计出单个或多个信号源的方位信息。由于激励信号和方向图之间存在傅里叶关系,DoA 估计也可以等效为谱估计问题。

DoA 估计:多重信号分类 MUSIC 算法(附 MATLAB 代码)

多重信号分类(Mutiple Signal Classification)算法,简称 MUSIC 算法,是一种常用的 DoA 估计方法。它的基本思想是将任意阵列输出数据的协方差矩阵进行特征分解,从而得到与信号分量相对应的信号子空间和与信号分量相正交的噪声子空间

信号模型

假设空间中存在 M M M 个不同方向的信号,入射到由 N N N 个天线单元构成的均匀直线阵上。第 i i i 个信号源的方向为 ϕ i \phi_i ϕi i = 1 , … , M i = 1,\dots,M i=1,,M),第 i i i 个信号源的信号为 α i ( t ) \alpha_i(t) αi(t)假设 M < N M < N M<N

令第 n n n 个天线单元的噪声 n n ( t ) n_n(t) nn(t)。在窄带远场条件下,第 n n n 个天线单元的输出信号 x n ( t ) x_n(t) xn(t) 可表示为
x n ( t ) = ∑ i = 1 M α i ( t ) e j k z n sin ⁡ ϕ i + n n ( t ) = ∑ i = 1 M α i ( t ) s n ( ϕ i ) + n n ( t ) \begin{aligned} x_n(t) &= \sum\limits_{i=1}^{M} \alpha_i(t) e^{jkz_n\sin\phi_i} + n_n(t) \\ &= \sum\limits_{i=1}^{M} \alpha_i(t) s_n(\phi_i) + n_n(t) \end{aligned} xn(t)=i=1Mαi(t)ejkznsinϕi+nn(t)=i=1Mαi(t)sn(ϕi)+nn(t)

N N N 个天线的输出信号表示成向量形式 x ( t ) \bf x (t) x(t),则上式可归纳为
x ( t ) = S α ( t ) + n ( t ) \mathbf{x}(t) = \mathbf{S} \mathbf{\alpha}(t) + \mathbf{n}(t) x(t)=Sα(t)+n(t)

其中, S \mathbf{S} S 为阵列的流型矩阵,矩阵规模为 N × M N\times M N×M,具体可表示为 M M M
个不同方向对应的阵列导向矢量
S = [ s ( ϕ 1 ) , s ( ϕ 2 ) , … , s ( ϕ M ) ] \mathbf{S} = [\mathbf{s}(\phi_1), \mathbf{s}(\phi_2), \dots, \mathbf{s}(\phi_M) ] S=[s(ϕ1),s(ϕ2),,s(ϕM)]

由于 M < N M < N M<N,流型矩阵 S \mathbf{S} S列满秩矩阵 R a n k ( S ) = M \mathrm{Rank}(\mathbf{S}) = M Rank(S)=M

MUSIC 算法思想

假设不同信号源的信号之间是相互正交的,噪声与信号之间是正交的,则阵列输出信号 x ( t ) \mathbf x(t) x(t) 的协方差矩阵
R = E [ x ( t ) x ( t ) H ] = E [ S α ( t ) α ( t ) H S H + n ( t ) n ( t ) H ] = S A S H + σ 2 I = R S + σ 2 I \begin{aligned} \mathbf R &= \mathrm E [\mathbf x(t) \mathbf x(t)^H] \\ &= \mathrm E [\mathbf{S} \mathbf{\alpha}(t) \mathbf{\alpha}(t)^H \mathbf{S}^H + \mathbf{n}(t)\mathbf{n}(t)^H ] \\ &= \mathbf{S} \mathbf A \mathbf{S}^H + \sigma^2 \mathbf I \\ &= \mathbf{R}_S + \sigma^2 \mathbf I \\ \end{aligned} R=E[x(t)x(t)H]=E[Sα(t)α(t)HSH+n(t)n(t)H]=SASH+σ2I=RS+σ2I

其中, A \mathbf A A 为不同信号源之间的协方差矩阵,由于不同信号源之间是相互正交的, A \mathbf A A正定对角矩阵
A = [ E [ ∣ α 1 ( t ) ∣ 2 ] … … … … E [ ∣ α 2 ( t ) ∣ 2 ] … … … … … … … … … E [ ∣ α M ( t ) ∣ 2 ] ] \mathbf A = \left[\begin{matrix} &\mathrm E [\mathbf |\alpha_1(t)|^2] &\dots &\dots &\dots \\ &\dots &\mathrm E [\mathbf |\alpha_2(t)|^2] &\dots &\dots \\ &\dots &\dots &\dots &\dots \\ &\dots &\dots &\dots &\mathrm E [\mathbf |\alpha_M(t)|^2] \\ \end{matrix}\right] A= E[α1(t)2]E[α2(t)2]E[αM(t)2]

由于信号协方差矩阵 R S \mathbf R_S RS 的规模为 N × N N\times N N×N秩为 M M M R S \mathbf R_S RS 存在 N − M N - M NM特征值为 0 的特征向量,令这种特征向量为 q m \mathbf q_m qm,则
R S q m = 0 \mathbf R_S \mathbf q_m = 0 RSqm=0
⇒ S A S H q m = 0 \Rightarrow \mathbf{S} \mathbf A \mathbf{S}^H \mathbf q_m = 0 SASHqm=0
⇒ q m H S A S H q m = 0 \Rightarrow \mathbf q_m^H \mathbf{S} \mathbf A \mathbf{S}^H \mathbf q_m = 0 qmHSASHqm=0
⇒ S H q m = 0 \Rightarrow \mathbf{S}^H \mathbf q_m = 0 SHqm=0

上述推论说明, R S \mathbf R_S RS 的特征值为 0 时对应的特征向量 q m \mathbf q_m qm 与信号源对应的 M M M 个导向矢量均正交

R S \mathbf R_S RS N − M N-M NM 个特征值为 0 时对应的特征向量构成矩阵 $\mathbf Q_n $,其规模为 N × ( N − M ) N \times (N-M) N×(NM),则
S H Q n = 0 \mathbf{S}^H \mathbf Q_n = 0 SHQn=0

MUSIC 算法的谱估计公式
P M U S I C ( ϕ ) = 1 ∣ ∣ Q n H s ( ϕ ) ∣ ∣ 2 P_{\mathrm{MUSIC}}(\phi) = \frac{1}{|| \mathbf Q_n^H \mathbf s(\phi) ||^2} PMUSIC(ϕ)=∣∣QnHs(ϕ)21

当上式中的 ϕ \phi ϕ 与信号源方向相同时,分母为零,此时 MUSIC 谱估计为无穷大。因此,MUSIC 谱估计的尖峰数目与信源数目相同,尖峰对应的方向即为信号源的方向

如何根据阵列输出信号 x \mathbf x x 计算 Q n \mathbf Q_n Qn​ ?

通过记录多组阵列输出信号快拍,可以计算出输出信号协方差矩阵的近似值
R = 1 K ∑ k = 1 K x k x k H \mathbf R = \frac{1}{K} \sum\limits_{k=1}^K \mathbf x_k \mathbf x_k^H R=K1k=1KxkxkH

那么,如何根据输出信号的协方差矩阵 R \mathbf R R 估计出信号协方差矩阵 R S \mathbf R_S RS 对应的特征值为 0 的特征向量矩阵 Q n \mathbf Q_n Qn 呢?

对于 R S \mathbf R_S RS 的任意特征向量 q m ∈ Q \mathbf q_m \in \mathbf Q qmQ ,有
R S q m = λ m q m \mathbf R_S \mathbf q_m = \lambda_m \mathbf q_m RSqm=λmqm
⇒ R q m = R S q m + σ 2 I q m = ( λ m + σ 2 ) q m \begin{aligned} \Rightarrow \mathbf R \mathbf q_m &= \mathbf R_S \mathbf q_m + \sigma^2 \mathbf I \mathbf q_m \\ &= (\lambda_m + \sigma^2)\mathbf q_m \end{aligned} Rqm=RSqm+σ2Iqm=(λm+σ2)qm

因此,信号协方差矩阵 R S \mathbf R_S RS 的特征值 λ m \lambda_m λm 对应的特征向量与输出信号协方差矩阵 R \mathbf R R 的特征值 λ m + σ 2 \lambda_m+\sigma^2 λm+σ2 对应的特征向量相同

因此, R \mathbf R R 的特征分解可表示为
R = Q ( Λ + σ 2 I ) Q H = Q [ λ 1 + σ 2 0 … 0 0 … 0 0 λ 2 + σ 2 … 0 0 … 0 … … … … … … … 0 0 … λ M + σ 2 0 … 0 0 0 … 0 σ 2 … 0 … … … … … … … 0 0 … 0 0 … σ 2 ] Q H \begin{aligned} \mathbf R &= \mathbf Q (\mathbf \Lambda + \sigma^2 \mathbf I) \mathbf Q^H \\ &= \mathbf Q \left[\begin{matrix} &\lambda_1 + \sigma^2 &0 &\dots &0 &0 &\dots &0 \\ &0 &\lambda_2 + \sigma^2 &\dots &0 &0 &\dots &0 \\ &\dots &\dots &\dots &\dots &\dots &\dots &\dots \\ &0 &0 &\dots &\lambda_M + \sigma^2 &0 &\dots &0 \\ &0 &0 &\dots &0 &\sigma^2 &\dots &0 \\ &\dots &\dots &\dots &\dots &\dots &\dots &\dots \\ &0 &0 &\dots &0 &0 &\dots &\sigma^2 \\ \end{matrix}\right] \mathbf Q^H \end{aligned} R=Q(Λ+σ2I)QH=Q λ1+σ200000λ2+σ200000λM+σ200000σ200000σ2 QH

上式表明,将输出信号矩阵 R \mathbf R R 进行特征分解,得到的 N − M N-M NM 个较小且相等的特征值对应的特征向量即可构成 Q n \mathbf Q_n Qn

MATLAB 仿真

DoA 估计:多重信号分类 MUSIC 算法(附 MATLAB 代码)

clc; clear; close all;

%% 参数设置
%%% 工作频率
c = 3e8;
freq = 10e9;
lambda = c/freq;    % 波长
k = 2*pi/lambda;    % 波数
%%% 阵列参数
N = 10;                 % 阵元数量
d = 0.5*lambda;         % 阵元间隔 
z = (0:d:(N-1)*d)';     % 阵元坐标分布
%%% 信号源参数
phi = [-10, -30, 60]'*pi/180;   % 来波方向
M = length(phi);                % 信号源数目
%%% 仿真参数
SNR = 10;             % 信噪比(dB)
K = 1000;     % 采样点数

%% 阵列接收信号仿真模拟
S = exp(1j*k*z*sin(phi'));          % 流型矩阵
Alpha = randn(M, K);         % 输入信号
X = S*Alpha;                        % 阵列接收信号
X1 = awgn(X, SNR, 'measured');      % 加载高斯白噪声

%% MUSIC 算法
%%% 阵列接收信号的协方差矩阵的特征分解
R = X1*X1'/K;    % 阵列接收信号的协方差矩阵
[EV, D] = eig(R);       % 特征值分解
EVA = diag(D);          % 提取特征值
[EVA, I] = sort(EVA, 'descend');   % 降序排序
Q = EV(:, I);           % 特征向量构成的矩阵
Q_n = Q(:, M+1:N);      % 噪声子空间
%%% 计算MUSIC谱估计函数
phi_list = linspace(-pi/2, pi/2, 200)';
S1 = exp(1j*k*z*sin(phi_list'));    % 不同方向对应的流型矢量构成矩阵
P_MUSIC = 1./sum(abs(Q_n'*S1).^2);     % MUSIC 谱估计公式
%%% 转换为dB
P_MUSIC = abs(P_MUSIC);
P_MUSIC_max = max(P_MUSIC);
P_MUSIC_dB = 10*log10(P_MUSIC/P_MUSIC_max);
%%% 提取信号源方向
[P_peaks, P_peaks_idx] = findpeaks(P_MUSIC_dB);     % 提取峰值
[P_peaks, I] = sort(P_peaks, 'descend');    % 峰值降序排序
P_peaks_idx = P_peaks_idx(I);
P_peaks = P_peaks(1:M);             % 提取前M个
P_peaks_idx = P_peaks_idx(1:M);
phi_e = phi_list(P_peaks_idx)*180/pi;   % 估计方向
disp('信号源估计方向为:');
disp(phi_e);
%%% 绘图
figure;
plot(phi_list*180/pi, P_MUSIC_dB, 'k', 'Linewidth', 2);
xlabel('\phi (deg)');
ylabel('Pseudo-spectrum (dB)');
axis([-90, 90, -40, 0]);
grid on;
hold on;
plot(phi_e, P_peaks, 'r.', 'MarkerSize', 25);
hold on;
for idx = 1:M
    text(phi_e(idx)+3, P_peaks(idx), sprintf('%0.1f°', phi_e(idx)));
end

参考文献

[1] 王永良. 空间谱估计理论与算法[M]. 清华大学出版社, 2004.
[2] 张小飞, 陈华伟, 仇小锋. 阵列信号处理及MATLAB实现[M]. 电子工业出版社, 2015.文章来源地址https://www.toymoban.com/news/detail-436088.html

到了这里,关于DoA 估计:多重信号分类 MUSIC 算法(附 MATLAB 代码)的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • DOA算法之DBF、CAPON、MUSIC、ROOT-MUSIC、ESPRIT、DML算法对比

    阵列信号处理算法应用领域涉及雷达、声纳、卫星通信等众多领域,其 主要目的就是对天线阵列接收到的信号进行处理,增强有用信号,抑制无用信号,以达到空域滤波的目的 ,最后提取回波信号中所包含的角度等信息。 DOA估计意思是波达角度估计,是指电磁波到达天线阵

    2024年02月06日
    浏览(130)
  • 【学习笔记】【DOA子空间算法】3 root-MUSIC 算法

      root-MUSIC 算法是 MUSIC 算法的一种多项式求根形式,其基本思想是 Pisarenko 分解。相比于 MUSIC 算法,root-MUSIC 算法无须谱峰搜索,降低了复杂度。   由前面可知,ULA 的方向矢量 a ( θ ) mathbf{a}(theta) a ( θ ) 如下: a ( θ ) = [ exp ⁡ ( − j 2 π d 1 sin ⁡ θ / λ ) ⋮ exp ⁡ ( −

    2024年02月12日
    浏览(36)
  • 【学习笔记】【DOA子空间算法】5 SS-MUSIC 算法

      在学习 SS-MUSIC 算法之前需要理解相干信号的概念。首先定义两个平稳信号 s i ( t ) s_i(t) s i ​ ( t ) 和 s k ( t ) s_k(t) s k ​ ( t ) 的相关系数 ρ i k rho_{ik} ρ ik ​ 如下: ρ i k = E [ s i ( t ) s k ∗ ( t ) ] E [ ∣ s i ( t ) ∣ 2 ] E [ ∣ s k ( t ) ∣ 2 ] begin{equation*} rho_{ik} = frac{mathrm{E}

    2024年02月12日
    浏览(59)
  • DOA估计算法——Capon算法

            在理解Capon算法之前,我们有必要先了解波束形成的基本思想以及原理到底是什么。这有助于我们更好的理解Capon算法的思想。 图 1  如图1展示了均匀阵列波束导向的示意图。图中 wm 表示加权值,波速形成(DBF)的基本思想就是将各阵元输出进行加权求和,在一定时

    2024年02月02日
    浏览(43)
  • DOA估计算法——迭代自适应算法(IAA)

            迭代自适应法 (Iterative Adaptive Approach,IAA)估计算法最早由美国的电气工程师和数学家Robert Schmidt和Roy A. Kuc在1986年的一篇论文\\\"Multiple Emitter Location and Signal Parameter Estimation\\\"中首次提出了这一算法, IAA DOA 估计算法是一种用于无线通信和雷达系统中估计信号到达方向的

    2024年01月21日
    浏览(36)
  • MUSIC算法相关原理知识(物理解读+数学推导+Matlab代码实现)

    部分来自于网络教程,如有侵权请联系本人删除  教程链接:MUSIC算法的直观解释:1,MUSIC算法的背景和基础知识_哔哩哔哩_bilibili  MUSIC算法的直观解释:2,我对于MUSIC算法的理解_哔哩哔哩_bilibili https://blog.csdn.net/zhangziju/article/details/100730081  一、MUSIC算法作用 MUSIC (Multiple

    2024年02月02日
    浏览(40)
  • 现代信号处理实验:MATLAB实现LD算法进行AR估计

    利用给定的一组样本数据估计一个平稳随机信号的功率谱密度称为功率谱估计,又称谱估计。谱估计的方法可以分成经典谱估计和现代谱估计。 经典谱估计又称为非参数化的谱估计,分为直接法和间接法。直接法是指直接计算样本数据的傅里叶变换,即获取频谱,然后计算频

    2024年02月03日
    浏览(55)
  • 稀疏DOA估计的经典算法——l1-SVD算法

    文献\\\"A Sparse Signal Reconstruction Perspective for Source Localization With Sensor Arrays\\\"提出了一种稀疏表示的DOA定位算法,它属于On-Grid类算法的范畴。其核心要点有二:其一是,通过了奇异值(SVD)分解,把以大量快拍数衡量的信号模型,转换成以信源数衡量的低维信号模型;其二是,以二阶

    2024年02月01日
    浏览(46)
  • 【MATLAB源码-第139期】基于matlab的OFDM信号识别与相关参数的估计,高阶累量/小波算法调制识别,循环谱估计,带宽估计,载波数目估计等等。

    在现代无线通信系统中,正交频分复用(OFDM)因其高效的频谱利用率、强大的抗多径衰落能力以及灵活的带宽分配等优势,成为了一种非常重要的调制技术。然而,随着无线通信网络的复杂性增加,对OFDM信号的识别与参数估计提出了更高的要求。这不仅是为了提高通信质量

    2024年02月19日
    浏览(44)
  • 关于克拉美罗下界(CRLB)-及不同DOA估计算法下的方差(性能)对比

        参数估计 在科研、工程乃至生活中都有广泛的应用。参数估计要解决的问题简单来说就是:基于一组观测数据,通过某种方法来获得我们想要的,与观测数据相关的一个或多个参数。     克拉美-罗界(Cramr-Rao Bound, CRB) 是 无偏估计 里我们常用的且十分重要的 一种对不同

    2024年04月13日
    浏览(119)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包