有没有房建设计的网站,wordpress 登录页,logo制作教程,设计工作室网站unitary MUSIC 算法 论文 A Unitary Transformation Method for Angle-of-Arrival Estimation 中提出了 unitary MUSIC 的算法#xff0c;直译就是酉 MUSIC 算法#xff0c;即酉变换 MUSIC 算法。该算法的目的是简化计算复杂度#xff0c;将传统 MUSIC 算法中的复数 SVD 和复…unitary MUSIC 算法 论文 A Unitary Transformation Method for Angle-of-Arrival Estimation 中提出了 unitary MUSIC 的算法直译就是酉 MUSIC 算法即酉变换 MUSIC 算法。该算法的目的是简化计算复杂度将传统 MUSIC 算法中的复数 SVD 和复数网格搜索计算转化为实数计算。在学习 unitary MUSIC 之前需要理解 Hermitian 矩阵及 Persymmetric 矩阵的概念及性质
Hermitian 矩阵指的是满足 A H A \mathbf{A}^H \mathbf{A} AHA 的矩阵 A \mathbf{A} APersymmetric 矩阵指的是满足 A J J A T \mathbf{A}\mathbf{J} \mathbf{J}\mathbf{A}^T AJJAT 的矩阵 A \mathbf{A} A其中 J \mathbf{J} J 其对角线从左下至右上很多地方又称 J \mathbf{J} J 为选择矩阵。假若矩阵 A \mathbf{A} A 既为 Hermitian 矩阵又为 Persymmetric 矩阵则满足 J A ∗ J A \mathbf{J}\mathbf{A}^*\mathbf{J}\mathbf{A} JA∗JA 其中 A ∗ \mathbf{A}^* A∗ 为 A \mathbf{A} A 的共轭。 在接下来的讨论中 I \mathbf{I} I 和 J \mathbf{J} J 分别用作表示单位矩阵和选择矩阵下文中将会出现这两种矩阵的运算例如 A I \mathbf{A}\mathbf{I} AI 或 J B \mathbf{J}\mathbf{B} JB设 A \mathbf{A} A 和 B \mathbf{B} B 均为方阵如果没有特别强调则说明 I \mathbf{I} I 和 J \mathbf{J} J 分别和 A \mathbf{A} A 和 B \mathbf{B} B 同维度。
算法原理 前面讨论的子空间算法中复协方差矩阵的特征值分解是至关重要的一步然而该步的计算量很高。为了降低计算量unitary MUSIC 算法考虑利用一个酉矩阵将原先的复协方差矩阵 R \mathbf{R} R 转换成实协方差矩阵同时传统算法中的复空间搜索向量 a ( θ ) \mathbf{a}(\theta) a(θ) 也用实向量来代替。 unitary MUSIC 算法的提出基于一个性质即若不相关的窄带远场信号源射入均匀线阵中其协方差矩阵不仅是 Hermitian且 Persymmetric。通常估计的协方差矩阵 R ≜ R ^ \mathbf{R}\triangleq \hat{\mathbf{R}} R≜R^ 仅仅只是 Hermitian 矩阵但不满足 Persymmetric 性质需要先获得一个满足 Persymmetric 性质的估计协方差矩阵 R ≜ 1 2 ( R ^ J R ^ ∗ J ) \begin{equation*} \mathbf{R}\triangleq \frac{1}{2}(\hat{\mathbf{R}} \mathbf{J}\hat{\mathbf{R}}^*\mathbf{J}) \end{equation*} R≜21(R^JR^∗J) 假设阵元数 M M M 为偶数unitary MUSIC 算法引入了一个酉矩阵 Q ∈ C M × M \mathbf{Q}\in\mathbb{C}^{M\times M} Q∈CM×M Q 1 2 [ I J j J − j I ] \mathbf{Q} \frac{1}{\sqrt{2}} \begin{bmatrix} \mathbf{I} \mathbf{J} \\ \mathrm{j}\mathbf{J} -\mathrm{j}\mathbf{I} \end{bmatrix} Q2 1[IjJJ−jI] 其中 I \mathbf{I} I 和 J \mathbf{J} J 分别为单位矩阵和选择矩阵且该两个矩阵维度均为 M 2 × M 2 \frac{M}{2}\times \frac{M}{2} 2M×2M。易得 Q \mathbf{Q} Q 为酉矩阵即 Q − 1 Q H \mathbf{Q}^{-1} \mathbf{Q}^H Q−1QH同时满足 Q ∗ J Q \mathbf{Q}^*\mathbf{J} \mathbf{Q} Q∗JQ 至此到了本算法的关键它在于证明由于 R \mathbf{R} R 是 Hermitian 且 Persymmetric 矩阵 Q R Q H \mathbf{Q}\mathbf{R}\mathbf{Q}^H QRQH 是实对称矩阵 因为 R \mathbf{R} R 为 Hermitian易得 Q R Q H \mathbf{Q}\mathbf{R}\mathbf{Q}^H QRQH 为 Hermitian因此只需要证明 Q R Q H \mathbf{Q}\mathbf{R}\mathbf{Q}^H QRQH 是实矩阵即证明 ( Q R Q H ) ∗ Q R Q H (\mathbf{Q}\mathbf{R}\mathbf{Q}^H)^* \mathbf{Q}\mathbf{R}\mathbf{Q}^H (QRQH)∗QRQH ( Q R Q H ) ∗ Q ∗ R ∗ Q T ( Q ∗ J ) ( J R ∗ J ) ( J Q T ) Q R Q H \begin{equation*} \begin{aligned} (\mathbf{Q}\mathbf{R}\mathbf{Q}^H)^* \\ \mathbf{Q}^*\mathbf{R}^*\mathbf{Q}^T \\ (\mathbf{Q}^*\mathbf{J})(\mathbf{J}\mathbf{R}^*\mathbf{J})(\mathbf{J}\mathbf{Q}^T)\\ \mathbf{Q}\mathbf{R}\mathbf{Q}^H \end{aligned} \end{equation*} (QRQH)∗Q∗R∗QT(Q∗J)(JR∗J)(JQT)QRQH 由此得证。 综上所述unitary MUSIC 算法引入酉矩阵 Q \mathbf{Q} Q 并令 R ≜ Q R Q H \mathbf{R}\triangleq\mathbf{Q}\mathbf{R}\mathbf{Q}^H R≜QRQH使得酉变换后的协方差矩阵变为实对称矩阵接着对其特征值分解即可进行后续的搜索步骤。而 ULA 的搜索方向矢量为 a ( θ ) [ 1 , e − j 2 π d sin θ / λ , ⋯ , e − j ( M − 1 ) 2 π d sin θ / λ ] T \mathbf{a}(\theta) \left[1, e^{-\mathrm{j}2\pi d\sin\theta/\lambda},\cdots,e^{-\mathrm{j}(M-1)2\pi d\sin\theta/\lambda}\right]^T a(θ)[1,e−j2πdsinθ/λ,⋯,e−j(M−1)2πdsinθ/λ]T 则 unitary MUSIC 算法的搜索方向矢量为 a ( θ ) ≜ Q a ( θ ) \mathbf{a}(\theta)\triangleq\mathbf{Q}\mathbf{a}(\theta) a(θ)≜Qa(θ)。 为了进一步降低算法计算复杂度Unitary MUSIC 算法考虑将搜索方向矢量也用实变量代替做法如下 a ( θ ) ≜ e j M − 1 2 2 π d sin θ / λ Q a ( θ ) Q [ e j M − 1 2 2 π d sin θ / λ , ⋯ , e j 1 2 2 π d sin θ / λ , e − j 1 2 2 π d sin θ / λ , ⋯ , e − j M − 1 2 2 π d sin θ / λ ] T \begin{aligned} \mathbf{a}(\theta) \triangleq e^{j\frac{M-1}{2}2\pi d\sin\theta/\lambda}\mathbf{Q}\mathbf{a}(\theta) \\ \mathbf{Q}\left[e^{\mathrm{j}\frac{M-1}{2}2\pi d\sin\theta/\lambda},\cdots, e^{\mathrm{j}\frac{1}{2}2\pi d\sin\theta/\lambda},e^{-\mathrm{j}\frac{1}{2}2\pi d\sin\theta/\lambda},\cdots,e^{-\mathrm{j}\frac{M-1}{2}2\pi d\sin\theta/\lambda}\right]^T \end{aligned} a(θ)≜ej2M−12πdsinθ/λQa(θ)Q[ej2M−12πdsinθ/λ,⋯,ej212πdsinθ/λ,e−j212πdsinθ/λ,⋯,e−j2M−12πdsinθ/λ]T 不难看出原方向矢量对应的阵列索引位置为 { 0 , 1 , ⋯ , M − 1 } \{0,1,\cdots,M-1\} {0,1,⋯,M−1}而更新后方向矢量对应的阵列索引位置为 { − M − 1 2 , − M − 3 2 , ⋯ , − 1 2 , 1 2 , ⋯ , M − 3 2 , M − 1 2 } \{-\frac{M-1}{2},-\frac{M-3}{2},\cdots,-\frac{1}{2},\frac{1}{2},\cdots,\frac{M-3}{2},\frac{M-1}{2}\} {−2M−1,−2M−3,⋯,−21,21,⋯,2M−3,2M−1}。此时 unitary MUSIC 算法的搜索方向矢量更新为 a ‾ ( θ ) e j M − 1 2 2 π d sin θ / λ Q a ( θ ) 2 [ cos ( M − 1 2 2 π d sin θ / λ ) ⋮ cos ( 1 2 2 π d sin θ / λ ) sin ( − 1 2 2 π d sin θ / λ ) ⋮ sin ( − M − 1 2 2 π d sin θ / λ ) ] \begin{aligned} \overline{\mathbf{a}}(\theta) e^{\mathrm{j}\frac{M-1}{2}2\pi d\sin\theta/\lambda}\mathbf{Q}\mathbf{a}(\theta) \\ \sqrt{2} \begin{bmatrix} \cos\left(\frac{M-1}{2}2\pi d\sin\theta/\lambda\right) \\ \vdots \\ \cos\left(\frac{1}{2}2\pi d\sin\theta/\lambda\right) \\ \sin\left(-\frac{1}{2}2\pi d\sin\theta/\lambda\right) \\ \vdots \\ \sin\left(-\frac{M-1}{2}2\pi d\sin\theta/\lambda\right) \end{bmatrix} \end{aligned} a(θ)ej2M−12πdsinθ/λQa(θ)2 cos(2M−12πdsinθ/λ)⋮cos(212πdsinθ/λ)sin(−212πdsinθ/λ)⋮sin(−2M−12πdsinθ/λ) 当 M M M 为奇数时酉矩阵 Q \mathbf{Q} Q 的形式为 Q 1 2 [ I O J O T 2 O T j J O − j I ] \mathbf{Q} \frac{1}{\sqrt{2}} \begin{bmatrix} \mathbf{I} \mathbf{O} \mathbf{J} \\ \mathbf{O}^T \sqrt{2} \mathbf{O}^T \\ \mathrm{j}\mathbf{J} \mathbf{O} -\mathrm{j}\mathbf{I} \end{bmatrix} Q2 1 IOTjJO2 OJOT−jI 其中 O \mathbf{O} O 为零矩阵。
进一步理解 在上一小节中我整理了论文对于 unitary MUSIC 算法的解释及证明其核心思想在于酉矩阵 Q \mathbf{Q} Q 的提出。在本小节中我将进一步对 Q \mathbf{Q} Q 的作用谈谈自己的理解。 在论文中unitary 算法一直强调 R \mathbf{R} R 的 Hermitian 及 Persymmetric 性质因为假若 R \mathbf{R} R 不满足这两个性质 Q R Q H \mathbf{Q}\mathbf{R}\mathbf{Q}^H QRQH 将不是实矩阵但是利用 R e ( Q R Q H ) \mathrm{Re}(\mathbf{Q}\mathbf{R}\mathbf{Q}^H) Re(QRQH) 来进行后续的估计其实是可以估计角度的。假设 M 4 M4 M4 和 φ 2 π d sin θ / λ \varphi 2\pi d \sin\theta/\lambda φ2πdsinθ/λ展开 Q a ( θ ) \mathbf{Q}\mathbf{a}(\theta) Qa(θ) Q a ( θ ) 2 [ 1 0 0 1 0 1 1 0 0 j − j 0 j 0 0 − j ] [ 1 e − j φ e − j 2 φ e − j 3 φ ] 2 [ 1 e − j 3 φ e − j φ e − j 2 φ j ( e − j φ − e − j 2 φ ) j ( 1 − e − j 3 φ ) ] 2 e − j 3 2 φ [ e j 3 2 φ e − j 3 2 φ e j 1 2 φ e − j 1 2 φ j ( e j 1 2 φ − e − j 1 2 φ ) j ( e j 3 2 φ − e − j 3 2 φ ) ] 2 e − j 3 2 φ [ cos ( 3 2 φ ) cos ( 1 2 φ ) sin ( − 1 2 φ ) sin ( − 3 2 φ ) ] e − j 3 2 φ a ‾ ( θ ) \begin{aligned} \mathbf{Q}\mathbf{a}(\theta) \sqrt{2} \begin{bmatrix} 1 0 0 1 \\ 0 1 1 0 \\ 0 \mathrm{j} -\mathrm{j} 0 \\ \mathrm{j} 0 0 -\mathrm{j} \end{bmatrix} \begin{bmatrix} 1\\ e^{-\mathrm{j}\varphi}\\ e^{-\mathrm{j}2\varphi}\\ e^{-\mathrm{j}3\varphi} \end{bmatrix} \\ \sqrt{2}\begin{bmatrix} 1 e^{-\mathrm{j}3\varphi}\\ e^{-\mathrm{j}\varphi} e^{-\mathrm{j}2\varphi}\\ \mathrm{j}(e^{-\mathrm{j}\varphi} - e^{-\mathrm{j}2\varphi})\\ \mathrm{j}(1 - e^{-\mathrm{j}3\varphi}) \end{bmatrix}\\ \sqrt{2}e^{-\mathrm{j}\frac{3}{2}\varphi}\begin{bmatrix} e^{\mathrm{j}\frac{3}{2}\varphi} e^{-\mathrm{j}\frac{3}{2}\varphi}\\ e^{\mathrm{j}\frac{1}{2}\varphi} e^{-\mathrm{j}\frac{1}{2}\varphi}\\ \mathrm{j}(e^{\mathrm{j}\frac{1}{2}\varphi} - e^{-\mathrm{j}\frac{1}{2}\varphi})\\ \mathrm{j}(e^{\mathrm{j}\frac{3}{2}\varphi} - e^{-\mathrm{j}\frac{3}{2}\varphi}) \end{bmatrix}\\ \sqrt{2}e^{-\mathrm{j}\frac{3}{2}\varphi}\begin{bmatrix} \cos(\frac{3}{2}\varphi) \\ \cos(\frac{1}{2}\varphi) \\ \sin(-\frac{1}{2}\varphi) \\ \sin(-\frac{3}{2}\varphi) \end{bmatrix}\\ e^{-\mathrm{j}\frac{3}{2}\varphi}\overline{\mathbf{a}}(\theta) \end{aligned} Qa(θ)2 100j01j001−j0100−j 1e−jφe−j2φe−j3φ 2 1e−j3φe−jφe−j2φj(e−jφ−e−j2φ)j(1−e−j3φ) 2 e−j23φ ej23φe−j23φej21φe−j21φj(ej21φ−e−j21φ)j(ej23φ−e−j23φ) 2 e−j23φ cos(23φ)cos(21φ)sin(−21φ)sin(−23φ) e−j23φa(θ) 至此可以得到 Q a ( θ ) e − j M − 1 2 φ a ‾ ( θ ) \mathbf{Q}\mathbf{a}(\theta)e^{-\mathrm{j}\frac{M-1}{2}\varphi}\overline{\mathbf{a}}(\theta) Qa(θ)e−j2M−1φa(θ)进一步我们可以得到 Q A e − j M − 1 2 φ A ‾ \mathbf{Q}\mathbf{A}e^{-\mathrm{j}\frac{M-1}{2}\varphi}\overline{\mathbf{A}} QAe−j2M−1φA 其中 A ‾ ∈ R M × K \overline{\mathbf{A}}\in\mathbb{R}^{M\times K} A∈RM×K 是由 K K K 个形如 a ‾ ( θ ) \overline{\mathbf{a}}(\theta) a(θ) 的实向量组成的矩阵。最后一步我们可以得到 Q R Q H ( e − j M − 1 2 φ A ‾ ) R s ( e j M − 1 2 φ A ‾ T ) A ‾ R s A ‾ T \begin{aligned} \mathbf{Q}\mathbf{R}\mathbf{Q}^H\left(e^{-\mathrm{j}\frac{M-1}{2}\varphi}\overline{\mathbf{A}}\right)\mathbf{R}_s\left(e^{\mathrm{j}\frac{M-1}{2}\varphi}\overline{\mathbf{A}}^T\right)\\ \overline{\mathbf{A}}\mathbf{R}_s\overline{\mathbf{A}}^T \end{aligned} QRQH(e−j2M−1φA)Rs(ej2M−1φAT)ARsAT 因此有 R e ( Q R Q H ) A ‾ R e ( R s ) A ‾ T \mathrm{Re}(\mathbf{Q}\mathbf{R}\mathbf{Q}^H) \overline{\mathbf{A}}\mathrm{Re}(\mathbf{R}_s)\overline{\mathbf{A}}^T Re(QRQH)ARe(Rs)AT不难看出即使 R \mathbf{R} R 不满足 Hermitian 及 Persymmetric 性质仍然不会破坏正交性。 总的来说 Q \mathbf{Q} Q 的作用就是使得方向矢量转为实向量如此便可以利用协方差矩阵的实部进行后续的计算。
算法步骤 unitary MUSIC 算法步骤如下输入为阵列接收矩阵 X \mathbf{X} X
计算协方差矩阵 R 1 T X X H \mathbf{R} \frac{1}{T} \mathbf{X}\mathbf{X}^H RT1XXH 和酉矩阵 Q \mathbf{Q} Q接着以 R ≜ Q R Q H \mathbf{R}\triangleq\mathbf{Q}\mathbf{R}\mathbf{Q}^H R≜QRQH 更新协方差矩阵。对 R \mathbf{R} R 进行特征值分解并对特征值进行排序然后取得 M − K M-K M−K 个较小特征值对应的特征向量来组成噪声子空间 U n \mathbf{U}_n Un。以下式遍历 θ ∈ [ − 9 0 ∘ , 9 0 ∘ ] \theta \in [-90^{\circ}, 90^{\circ}] θ∈[−90∘,90∘] P ( θ ) 1 a ‾ H ( θ ) U n U n T a ‾ ( θ ) \begin{equation*} P(\theta) \frac{1}{\overline{\mathbf{a}}^H(\theta)\mathbf{U}_n\mathbf{U}_n^T\overline{\mathbf{a}}(\theta)} \end{equation*} P(θ)aH(θ)UnUnTa(θ)1 此时得到一组 P ( θ ) P(\theta) P(θ) K K K 个最大值对应的 θ \theta θ 就是需要返回的结果。
代码实现
clear; close all; clc;%% Parameters
lambda 3e8/1e9; % wavelength, c/f
d lambda/4; % distance between sensors
theta [10,20]; % true DoAs, 1 times K vector
theta sort(theta);
M 16; % # of sensors
T 500; % # of snapshots
K length(theta); % # of signals
noise_flag 1;
SNR 0; % signal-to-noise ratio
grid 0.1; % search grid%% Signals
R generateSignal(M,K,T,theta,lambda,d,noise_flag,SNR);%% DoA
% unitary-MUSIC
[theta_unitary_music,P_unitary_music] unitaryMUSIC(R,M,K,lambda,d,grid);%% plot
figure;
hold on;
ang_list -90:grid:90;
plot(ang_list, P_unitary_music);
hold off;function [R,X,A,S] generateSignal(M,K,T,theta,lambda,d,noise_flag,SNR)S exp(1j*2*pi*randn(K,T)); % signal matrixA exp(-1j*(0:M-1)*2*pi*d/lambda*sind(theta)); % steering vector matrixN noise_flag.*sqrt(10.^(-SNR/10))*(1/sqrt(2))*(randn(M,T)1j*randn(M,T)); % noise matrixX A*SN; % received matrixR X*X/T; % covariance matrix
endfunction [theta,P] unitaryMUSIC(R,M,K,lambda,d,grid)M_half floor(M/2);O zeros(M_half,1);I eye(M_half);J fliplr(I);if mod(M,2) 0Q [I,J;1j*J,-1j*I]./sqrt(2);elseQ [I,O,J;O,sqrt(2),O;1j*J,O,-1j*I]./sqrt(2);endR real(Q*R*Q);[U,~] svd(R);Un U(:, K1:M);a_list exp(-1j*(0:M-1).*2*pi*d/lambda*sind(-90:grid:90));a_list a_list.*exp(1j*(M-1)*ones(M,1)*pi*d/lambda*sind(-90:grid:90));a_list real(Q*a_list);P arrayfun((i) 1/norm(Un*a_list(:,i),fro),1:size(a_list,2)); % spectral spectrum grid searchP 10*log10(P./max(P));[~, idx] findpeaks(P,NPeaks,K,SortStr,descend); % find K peakstheta sort((idx-1)*grid-90);
end参考内容
Huarng K C, Yeh C C. A unitary transformation method for angle-of-arrival estimation[J]. IEEE Transactions on Signal Processing, 1991, 39(4): 975-977.【wikipedia】Persymmetric matrix【wikipedia】Hermitian matrix
