多重信號分類MUSIC演算法

1. MUSIC演算法原理

MUSIC演算法，叫做多信號分類演算法 (Multiple Signal Classification)，是一種基于特征結構的高解析度DOA演算法，該演算法利用了信號子空間和噪聲子空間正交性的特點，構造噪聲空間然后通過譜峰搜索來檢測信號的波達方向，需要注意的是，該演算法有一個前提，即各個入射信號之間互不相關，這樣才能保證入射信號的協方差矩陣是滿秩的，

考慮以下這個線陣模型：

令 X ( t ) X(t) X(t) 是在第 t t t 個快拍觀察到的資料向量，在陣列信號處理和空間譜估計中, X ( t ) = \boldsymbol{X}(t)= X(t)= [ X 1 ( t ) , X 2 ( t ) , ? ? , X N ( t ) ] T \left[X_{1}(t), X_{2}(t), \cdots, X_{N}(t)\right]^{\mathrm{T}} [X1?(t),X2?(t),?,XN?(t)]T 由 N N N 個陣元 (天線或傳感器) 的觀測資料組成，在時域譜估計中，向量 x ( t ) = [ x ( t ) , x ( t ? 1 ) , ? ? , x ( t ? n ? 1 ) ] T \boldsymbol{x}(t)=[x(t), x(t-1), \cdots, x(t-n-1)]^{\mathrm{T}} x(t)=[x(t),x(t?1),?,x(t?n?1)]T 由連續的 n n n 個觀察資料樣本組成，
假定資料向量 X ( t ) \boldsymbol{X}(t) X(t) 是 r r r 個窄帶信號入射到 N N N 個陣元組成的陣列的觀察資料向量或者是 r r r 個不相干的復諧波的疊加，即
X ( t ) = ∑ i = 1 r s i ( t ) a ( ω i ) + v ( t ) = A s ( t ) + n ( t ) \boldsymbol{X}(t)=\sum_{i=1}^{r} s_{i}(t) \boldsymbol{a}\left(\omega_{i}\right)+\boldsymbol{v}(t)=\boldsymbol{A} \boldsymbol{s}(t)+\boldsymbol{n}(t) X(t)=i=1∑r?si?(t)a(ωi?)+v(t)=As(t)+n(t)式中， A = [ a ( ω 1 ) , ? ? , a ( ω r ) ] \boldsymbol{A}=\left[\boldsymbol{a}\left(\omega_{1}\right), \cdots, \boldsymbol{a}\left(\omega_{r}\right)\right] A=[a(ω1?),?,a(ωr?)] 為 N × r N \times r N×r 陣列流型矩陣， a ( ω i ) = [ 1 , e j ω i , … , e j ( N ? 1 ) ω i ] T \boldsymbol{a}\left(\omega_{i}\right)=\left[1, \mathrm{e}^{\mathrm{j} \omega_{i}}, \ldots, \mathrm{e}^{\mathrm{j}(N-1) \omega_{i}}\right]^{\mathrm{T}} a(ωi?)=[1,ejωi?,…,ej(N?1)ωi?]T 為方向向量; s ( t ) = [ s 1 ( t ) , ? ? , s r ( t ) ] T \boldsymbol{s}(t)=\left[s_{1}(t), \cdots, s_{r}(t)\right]^{\mathrm{T}} s(t)=[s1?(t),?,sr?(t)]T 為隨機信號向量，其均值為零向量，協方差矩陣為 R s = E { s ( t ) s H ( t ) } ; \boldsymbol{R}_{\boldsymbol{s}}=\mathrm{E}\left\{\boldsymbol{s}(t) \boldsymbol{s}^{\mathrm{H}}(t)\right\} ; Rs?=E{s(t)sH(t)}; 而 v ( t ) = [ v 1 ( t ) , ? ? , v N ( t ) ] T \boldsymbol{v}(t)=\left[v_{1}(t), \cdots, v_{N}(t)\right]^{\mathrm{T}} v(t)=[v1?(t),?,vN?(t)]T 為加性噪聲向量，其各個分量為高斯白噪聲，它們具有零均值和相同的方差 σ 2 \sigma^{2} σ2 ，在諧波恢復中，引數 ω i \omega_{i} ωi? 為復諧波的頻率 ; 在陣列信號處理中, ω i \omega_{i} ωi? 是一空間引數
ω i = 2 π d λ sin ? θ i \omega_{i}=2 \pi \frac{d}{\lambda} \sin \theta_{i} ωi?=2πλd?sinθi?式中, d d d 為相鄰兩個陣元之間的距離 (假定陣元等間距排列成一直線)， λ \lambda λ 為波長, 且 θ i \theta_{i} θi? 表 示第 i i i 個窄帶信號達到陣元的入射方向，簡稱波達方向，
MUSIC演算法要解決的是根據 N N N 個快拍的觀測資料向量 X ( t ) ( t = 1 , 2 , ? ? , N ) \boldsymbol{X}(t)(t=1,2, \cdots, N) X(t)(t=1,2,?,N) 估計 r r r 個引數 ω i ° \omega_{i \circ} ωi°? 這相當于對 r r r 個混合信號進行分類，簡稱多重信號分類，
假定噪聲向量 v ( t ) \boldsymbol{v}(t) v(t) 與信號向量 s ( t ) \boldsymbol{s}(t) s(t) 統計不相關，并令觀測資料向量的協方差矩陣 R x x = E { x ( t ) x H ( t ) } \boldsymbol{R}_{x x}=\mathrm{E}\left\{\boldsymbol{x}(t) \boldsymbol{x}^{\mathrm{H}}(t)\right\} Rxx?=E{x(t)xH(t)} 的特征值分解為
R x x = A R s s A H + σ 2 I = U Σ U H = [ U S , U N ] [ Σ S O O Σ N ] [ U S H U N H ] \boldsymbol{R}_{x x}=\boldsymbol{A R}_{s s} \boldsymbol{A}^{\mathrm{H}}+\sigma^{2} \boldsymbol{I}=\boldsymbol{U} \boldsymbol{\Sigma} \boldsymbol{U}^{\mathrm{H}}=[\boldsymbol{U_S}, \boldsymbol{U_N}]\left[\begin{array}{cc} \boldsymbol{\Sigma_S} & \boldsymbol{O} \\ \boldsymbol{O} & \boldsymbol{\Sigma_N} \end{array}\right]\left[\begin{array}{l} \boldsymbol{U}^{\mathrm{H}}_S \\ \boldsymbol{U}^{\mathrm{H}}_N \end{array}\right] Rxx?=ARss?AH+σ2I=UΣUH=[US?,UN?][ΣS?O?OΣN??][USH?UNH??]式中, R s s = E { s ( t ) s H ( t ) } , \boldsymbol{R}_{s s}=\mathrm{E}\left\{\boldsymbol{s}(t) \boldsymbol{s}^{\mathrm{H}}(t)\right\}, Rss?=E{s(t)sH(t)}, 且 Σ S \boldsymbol{\Sigma_S} ΣS? 包含了 r r r 個大特征值, 它們比 σ 2 \sigma^{2} σ2 明顯大很多， Σ N = σ 2 I n ? r \boldsymbol{\Sigma_N}=\sigma^{2} I_{n-r} ΣN?=σ2In?r?，
因此有：
R x x U N = [ U S , U N ] [ Σ S O O Σ N ] [ U S H U N H ] U N = [ U S , U N ] [ Σ S O O Σ N ] [ O I ] = σ 2 U N \boldsymbol{R}_{x x} \boldsymbol{U_N}=[\boldsymbol{U_S}, \boldsymbol{U_N}]\left[\begin{array}{cc} \boldsymbol{\Sigma_S} & \boldsymbol{O} \\ \boldsymbol{O} & \boldsymbol{\Sigma_N} \end{array}\right]\left[\begin{array}{l} \boldsymbol{U}^{\mathrm{H}}_S \\ \boldsymbol{U}^{\mathrm{H}}_N \end{array}\right]\boldsymbol{U_N}=[\boldsymbol{U_S}, \boldsymbol{U_N}]\left[\begin{array}{cc}\boldsymbol{\Sigma_S} & \boldsymbol{O} \\ \boldsymbol{O} & \boldsymbol{\Sigma_N}\end{array}\right]\left[\begin{array}{l}\boldsymbol{O} \\ \boldsymbol{I}\end{array}\right]=\sigma^{2} \boldsymbol{U_N} Rxx?UN?=[US?,UN?][ΣS?O?OΣN??][USH?UNH??]UN?=[US?,UN?][ΣS?O?OΣN??][OI?]=σ2UN?
又由 R x x = A R s s A H + σ 2 I \boldsymbol{R}_{x x}=\boldsymbol{A R}_{s s} \boldsymbol{A}^{\mathrm{H}}+\sigma^{2} \boldsymbol{I} Rxx?=ARss?AH+σ2I 有 R x x U N = A R s s A H U N + σ 2 U N , \boldsymbol{R}_{x x} \boldsymbol{U_N}=\boldsymbol{A R}_{s s} \boldsymbol{A}^{\mathrm{H}} \boldsymbol{U_N}+\sigma^{2} \boldsymbol{U_N}, Rxx?UN?=ARss?AHUN?+σ2UN?,聯立上式可以得到
A R s s A H U N = O \boldsymbol{A R}_{s s} \boldsymbol{A}^{\mathrm{H}} \boldsymbol{U_N}=\boldsymbol{O} ARss?AHUN?=O進而有
U N H A R s s A H U N = O \boldsymbol{U}^{\mathrm{H}}_N \boldsymbol{A R}_{s s} \boldsymbol{A}^{\mathrm{H}} \boldsymbol{U_N}=\boldsymbol{O} UNH?ARss?AHUN?=O
眾所周知, Q \boldsymbol{Q} Q 非奇異時 t H Q t = 0 , \boldsymbol{t}^{\mathrm{H}} \boldsymbol{Q} \boldsymbol{t}=0, tHQt=0, 當且僅當 t = 0 , \boldsymbol{t}=\mathbf{0}, t=0, 故上式成立的充分必要條件是
A H U N = O \boldsymbol{A}^{\mathrm{H}} \boldsymbol{U_N}=\boldsymbol{O} AHUN?=O因為 R s s = E { s ( t ) s H ( t ) } R_{s s}=\mathrm{E}\left\{s(t) s^{\mathrm{H}}(t)\right\} Rss?=E{s(t)sH(t)} 非奇異，將 A = [ a ( ω 1 ) , ? ? , a ( ω p ) ] \boldsymbol{A}=\left[\boldsymbol{a}\left(\omega_{1}\right), \cdots, \boldsymbol{a}\left(\omega_{p}\right)\right] A=[a(ω1?),?,a(ωp?)] 代入上式即有
a H ( ω ) G = 0 T , ω = ω 1 , ω 2 , ? ? , ω p \boldsymbol{a}^{\mathrm{H}}(\omega) \boldsymbol{G}=\mathbf{0}^{\mathrm{T}}, \quad \omega=\omega_{1}, \omega_{2}, \cdots, \omega_{p} aH(ω)G=0T,ω=ω1?,ω2?,?,ωp?顯然, 當 ω ≠ ω 1 , ω 2 , ? ? , ω p \omega \neq \omega_{1}, \omega_{2}, \cdots, \omega_{p} ω?=ω1?,ω2?,?,ωp? 時, a H ( ω ) G ≠ 0 T \boldsymbol{a}^{\mathrm{H}}(\omega) \boldsymbol{G} \neq \mathbf{0}^{\mathrm{T}} aH(ω)G?=0T ，
將上式改寫成標量形式, 可以定義一種類似于功率譜的函式
P ( ω ) = 1 a H ( ω ) U N U N H a ( ω ) \boldsymbol P(\omega)=\frac{1}{\boldsymbol a^{\mathrm{H}}(\omega) \boldsymbol U_N \boldsymbol U^{\mathrm{H}}_N \boldsymbol a(\omega)} P(ω)=aH(ω)UN?UNH?a(ω)1?由于噪聲的影響， a H ( ω ) U N U N H a ( ω ) ≠ 0 \boldsymbol{a}^{\mathrm{H}}(\omega) \boldsymbol{U_N} \boldsymbol{U}^{\mathrm{H}}_N \boldsymbol{a}(\omega)\neq0 aH(ω)UN?UNH?a(ω)?=0，而是一個很小的數，因此此時的 P ( ω ) \boldsymbol P(\omega) P(ω)是一個極大值，對上式取峰值的 r r r 個 ω \omega ω 值 ω 1 , ω 2 , ? ? , ω r \omega_{1}, \omega_{2}, \cdots, \omega_{r} ω1?,ω2?,?,ωr? 給出 r r r 個信號的波達方向 θ 1 , θ 2 , ? ? , θ r \theta_{1}, \theta_{2}, \cdots, \theta_{r} θ1?,θ2?,?,θr?，
同時在實際應用中，接收陣接收到的資料有限長，我們可以通過這些接收資料估計入射信號的協方差矩陣，達到極大似然估計的目的，假設有M個快拍snapshot，則有
R ^ x x = 1 M ∑ i = 1 M x i x i H \hat{R}_{xx}=\frac{1}{M} \sum_{i=1}^{M} x_{i} x_{i}^{H} R^xx?=M1?i=1∑M?xi?xiH?然后對協方差矩陣 R ^ x x \hat{R}_{xx} R^xx?進行特征值分解：
R ^ x x = U Σ U H \begin{aligned} \hat{R}_{xx} &={U} \Sigma{U}^{H} \end{aligned} R^xx??=UΣUH?特征值分解后按照特征值的大小(從大到小)對特征向量進行重排，接下來就是要確定噪聲子空間，假設排序后的特征值 σ 1 2 ? σ 2 2 ? ? ? σ r 2 ? σ r + 1 2 ? ? ? σ N 2 \sigma_{1}^{2} \geqslant \sigma_{2}^{2} \geqslant \cdots \geqslant \sigma_{r}^{2} \geqslant \sigma_{r+1}^{2}\geqslant \cdots \geqslant \sigma_{N}^{2} σ12??σ22????σr2??σr+12????σN2?，接收信號總能量定義為 P = σ 1 2 + σ 2 2 + ? + σ N 2 P=\sigma_{1}^{2}+\sigma_{2}^{2}+\cdots+\sigma_{N}^{2} P=σ12?+σ22?+?+σN2?根據能量準則，如果J是滿足下式的最小整數： σ 1 2 + σ 2 2 + ? + σ J 2 P ? 0.95 或 0.9 \frac{\sigma_{1}^{2}+\sigma_{2}^{2}+\cdots+\sigma_{J}^{2}}{P}\geqslant0.95或0.9 Pσ12?+σ22?+?+σJ2???0.95或0.9那么 { σ J + 1 2 , σ J + 2 2 ? ? , σ N 2 } \{\sigma_{J+1}^{2},\sigma_{J+2}^{2}\cdots,\sigma_{N}^{2}\} {σJ+12?,σJ+22??,σN2?}對應的特征向量組成的空間 U N ^ \hat{U_N} UN?^?即為噪聲子空間，

然后就可以用 P ( ω ) = 1 a H ( ω ) U ^ N U ^ N H a ( ω ) \boldsymbol P(\omega)=\frac{1}{\boldsymbol a^{\mathrm{H}}(\omega) \boldsymbol {\hat U_N} \boldsymbol {\hat U^{\mathrm{H}}_N} \boldsymbol a(\omega)} P(ω)=aH(ω)U^N?U^NH?a(ω)1?進行譜峰搜索，尋找r個入射信號的來波方向，

2. MATLAB仿真代碼

% Code For Music Algorithm
% Author:癢羊羊
% Date：2020/10/28

clc; clear all; close all;
%% -------------------------initialization-------------------------
f = 500;                                        % frequency
c = 1500;                                       % speed sound
lambda = c/f;                                   % wavelength
d = lambda/2;                                   % array element spacing
M = 10;                                         % number of array elements
N = 100;                                        % number of snapshot
K = 6;                                          % number of sources
doa_phi = [-30, 0, 20, 40, 60, 75];             % direction of arrivals

%% generate signal
dd = (0:M-1)'*d;                                % distance between array elements and reference element
A = exp(-1i*2*pi*dd*sind(doa_phi)/lambda);      % manifold array, M*K
S = sqrt(2)\(randn(K,N)+1i*randn(K,N));         % array of random signal, K*N
X = A*S;                                        % received data without noise, M*N
X = awgn(X,10,'measured');                      % received data with SNR 10dB

%% calculate the covariance matrix of received data and do eigenvalue decomposition
Rxx = X*X'/N;                                   % covariance matrix
[U,V] = eig(Rxx);                               % eigenvalue decomposition
V = diag(V);                                    % vectorize eigenvalue matrix
[V,idx] = sort(V,'descend');                    % sort the eigenvalues in descending order
U = U(:,idx);                                   % reset the eigenvector
P = sum(V);                                     % power of received data
P_cum = cumsum(V);                              % cumsum of V

%% define the noise space
J = find(P_cum/P>=0.95);                        % or the coefficient is 0.9
J = J(1);                                       % number of principal component
Un = U(:,J+1:end);

%% music for doa; seek the peek
theta = -90:0.1:90;                             % steer theta
doa_a = exp(-1i*2*pi*dd*sind(theta)/lambda);    % manifold array for seeking peak
music = abs(diag(1./(doa_a'*(Un*Un')*doa_a)));  % the result of each theta
music = 10*log10(music/max(music));             % normalize the result and convert it to dB

%% plot
figure;
plot(theta, music, 'linewidth', 2);
title('Music Algorithm For Doa', 'fontsize', 16);
xlabel('Theta(°)', 'fontsize', 16);
ylabel('Spatial Spectrum(dB)', 'fontsize', 16);
grid on;

運行結果：

可以看到，入射信號互不相關時，傳統MUSIC演算法能夠以高解析度檢測出6個源的大致波達方向，分別是-29.7°，0°，19.8°，39.8°，60.4°，74.7°，但仍存在估計精度的問題，有很多改進的MUSIC演算法可以改善，
需要注意，陣元數目為M的半波長均勻線陣的自由度為DOF為M-1，表示該線陣最多可以分辨源的數目為M-1，同時若存在相干源，MUSIC演算法的效果就會不理想，基于相干源的MUSIC演算法會在后續文章詳寫，

部分MUSIC演算法原理參考張賢達先生的《矩陣分析與應用》，
歡迎轉載，表明出處，