root MUSIC 算法补充说明
root MUSIC 算法补充说明
- 多项式求根
- root MUSIC 算法原理
- 如何从 2 M − 2 2M-2 2M−2 个根中确定 K K K 个根
- 从复数域上观察
2
M
−
2
2M-2
2M−2 个根的分布
这篇笔记是上一篇关于 root MUSIC 笔记的补充。
多项式求根
要理解 root MUSIC 算法，需要理解多项式求根的相关知识。给定多项式 P ( x ) P(x) P(x)：
P ( x ) = a 0 + a 1 x + ⋯ + a n x n P(x) = a_0 + a_1 x + \cdots + a_n x^n P(x)=a0+a1x+⋯+anxn
容易看出 P ( x ) P(x) P(x) 中只有一个未知数 x x x，且未知数的最高次数为 n n n，因此称 P ( x ) P(x) P(x) 为一元 n n n 次多项式，同时系数 { a i ∈ C : i = 0 , ⋯ , n } \{a_i\in\mathbb{C}:i = 0,\cdots, n\} {ai∈C:i=0,⋯,n}。而多项式求根就是求得一元 n n n 次方程式 P ( x ) = 0 P(x)=0 P(x)=0 的解，这个解被称作根或者零点。
在进行后续的讨论前，还需要清楚，根据代数基本定理， n n n 次复系数多项式方程在复数域内有且只有 n n n 个根（这里的重根按重数计算）。
root MUSIC 算法原理
root MUSIC 算法是 MUSIC 算法的一种多项式求根形式。回忆一下，传统 MUSIC 算法利用了噪声子空间矩阵 U n \mathbf{U}_n Un 和搜索方向矢量 a ( θ ) \mathbf{a}(\theta) a(θ) 来构造空间谱：
P ( θ ) = 1 a H ( θ ) U n U n H a ( θ ) a ( θ ) = [ 1 , e − j 2 π d sin θ / λ , ⋯ , e − j 2 π ( M − 1 ) d sin θ / λ ] T P(\theta) = \frac{1}{\mathbf{a}^H(\theta)\mathbf{U}_n\mathbf{U}^H_n\mathbf{a}(\theta)} \\ \mathbf{a}(\theta) = \left[1,e^{-\mathrm{j}2\pi d \sin \theta/\lambda},\cdots,e^{-\mathrm{j}2\pi(M-1) d \sin \theta/\lambda}\right]^T P(θ)=aH(θ)UnUnHa(θ)1a(θ)=[1,e−j2πdsinθ/λ,⋯,e−j2π(M−1)dsinθ/λ]T
在 { θ = θ k : k = 1 , ⋯ , K } \{\theta = \theta_k:k = 1,\cdots,K\} {θ=θk:k=1,⋯,K} 时 P ( θ ) P(\theta) P(θ) 将产生峰值，换句话说此时 P − 1 ( θ ) = 0 P^{-1}(\theta)=0 P−1(θ)=0。
在接下来的讨论中，我们令 P − 1 ( θ ) = a H ( θ ) G a ( θ ) P^{-1}(\theta) = \mathbf{a}^H(\theta)\mathbf{G}\mathbf{a}(\theta) P−1(θ)=aH(θ)Ga(θ)，此时我们可以知道，MUSIC 算法满足 G ≜ U n U n H \mathbf{G} \triangleq \mathbf{U}_n\mathbf{U}^H_n G≜UnUnH，而 Capon 算法满足 G ≜ R − 1 \mathbf{G} \triangleq \mathbf{R}^{-1} G≜R−1。需要注意的是无论是 MUSIC 算法还是 Capon 算法， G \mathbf{G} G 均是 Hermitian 矩阵。
令 ω = − 2 π d sin θ / λ \omega = -2\pi d \sin\theta/\lambda ω=−2πdsinθ/λ 以及 z = e j ω z = e^{\mathrm{j}\omega} z=ejω，我们将会得到：
a ( z ) = [ 1 , z , z 2 , ⋯ , z M − 1 ] T = a ( θ ) P − 1 ( z ) = a H ( z ) G a ( z ) = P − 1 ( θ ) \begin{aligned} \mathbf{a}(z) &= [1,z,z^{2},\cdots,z^{M-1}]^T = \mathbf{a}(\theta) \\ P^{-1}(z) &= \mathbf{a}^H(z)\mathbf{G}\mathbf{a}(z) = P^{-1}(\theta) \end{aligned} a(z)P−1(z)=[1,z,z2,⋯,zM−1]T=a(θ)=aH(z)Ga(z)=P−1(θ)
接下来我们展开 P − 1 ( z ) P^{-1}(z) P−1(z)：
P − 1 ( z ) = a H ( z ) G a ( z ) = [ 1 , z ∗ , ( z ∗ ) 2 , ⋯ , ( z ∗ ) M − 1 ] G [ 1 , z , z 2 , ⋯ , z M − 1 ] T = [ 1 , z − 1 , z − 2 , ⋯ , z − M + 1 ] G [ 1 , z , z 2 , ⋯ , z M − 1 ] T = ∑ m = 0 M − 1 ∑ n = 0 M − 1 z − m G [ m , n ] z n = ∑ m = 0 M − 1 ∑ n = 0 M − 1 z n − m G [ m , n ] = ∑ p = − M + 1 M − 1 a p z − p \begin{aligned} P^{-1}(z) &= \mathbf{a}^H(z)\mathbf{G}\mathbf{a}(z) \\ &= [1,z^{*},(z^{*})^2,\cdots,(z^*)^{M-1}] \mathbf{G} [1,z,z^{2},\cdots,z^{M-1}]^T \\ &= [1,z^{-1},z^{-2},\cdots,z^{-M+1}] \mathbf{G} [1,z,z^{2},\cdots,z^{M-1}]^T \\ &= \sum_{m = 0}^{M-1} \sum_{n=0}^{M-1} z^{-m} \mathbf{G}_{[m,n]} z^{n} \\ &= \sum_{m = 0}^{M-1} \sum_{n=0}^{M-1} z^{n-m} \mathbf{G}_{[m,n]} \\ &=\sum_{p=-M+1}^{M-1}a_p z^{-p} \end{aligned} P−1(z)=aH(z)Ga(z)=[1,z∗,(z∗)2,⋯,(z∗)M−1]G[1,z,z2,⋯,zM−1]T=[1,z−1,z−2,⋯,z−M+1]G[1,z,z2,⋯,zM−1]T=m=0∑M−1n=0∑M−1z−mG[m,n]zn=m=0∑M−1n=0∑M−1zn−mG[m,n]=p=−M+1∑M−1apz−p
其中 G [ m , n ] \mathbf{G}_{[m,n]} G[m,n] 表示矩阵 G \mathbf{G} G 的第 m m m 行第 n n n 列元素， a p a_p ap 表示矩阵 G \mathbf{G} G 的第 p p p 条对角线的求和：
a p ≜ ∑ m − n = p G [ m , n ] a_p \triangleq \sum_{m-n = p} \mathbf{G}_{[m,n]} ap≜m−n=p∑G[m,n]
到这里我们已经可以看出，传统 MUSIC 算法对 P ( θ ) P(\theta) P(θ) 求峰值，其实等价于对 P − 1 ( z ) P^{-1}(z) P−1(z) 求根，为了方便大家的理解，我们令 M = 3 M=3 M=3，此时会得到一条简单的式子：
P − 1 ( z ) = a 2 z − 2 + a 1 z − 1 + a 0 z 0 + a − 1 z 1 + a − 2 z 2 P^{-1}(z) = a_{2}z^{-2}+a_{1}z^{-1} + a_{0}z^{0} + a_{-1}z^{1} + a_{-2}z^{2} P−1(z)=a2z−2+a1z−1+a0z0+a−1z1+a−2z2
可以看出，其实 P − 1 ( θ ) P^{-1}(\theta) P−1(θ) 是一个 2 M − 1 = 5 2M-1 = 5 2M−1=5 项的多项式，但还存在一个问题， P − 1 ( θ ) P^{-1}(\theta) P−1(θ) 中存在负整数次数，我们令 P − 1 ( z ) z M − 1 P^{-1}(z)z^{M-1} P−1(z)zM−1 将负整数次数消除即可，操作前后，求根的结果是一样的，因此我们可以说 P − 1 ( z ) z M − 1 P^{-1}(z)z^{M-1} P−1(z)zM−1 是一个一元的 2 M − 1 2M-1 2M−1 项的 2 M − 2 2M-2 2M−2 次的多项式。更进一步地，我们可以说求解 P − 1 ( z ) z M − 1 = 0 P^{-1}(z)z^{M-1}=0 P−1(z)zM−1=0 将会得到 2 M − 2 2M-2 2M−2 个根，从已知条件我们知道，其中 K K K 个根必定是 e j ω k e^{\mathrm{j}\omega_k} ejωk（ e j ω k e^{\mathrm{j}\omega_k} ejωk 的幅值是 1 1 1，因此该 K K K 点在单位圆上），在这里 ω k = − 2 π d sin θ k / λ \omega_k = -2\pi d \sin\theta_k/\lambda ωk=−2πdsinθk/λ。
总结一下，MUSIC 算法的谱峰搜索操作等价于对方程式 P − 1 ( z ) z M − 1 = 0 P^{-1}(z)z^{M-1}=0 P−1(z)zM−1=0 求根，root MUSIC 算法所做的，就是利用 G \mathbf{G} G 的多条对角线求和得到对应的多项式系数，从而求解得 2 M − 2 2M-2 2M−2 个根，接着筛选得到合适的 K K K 个根 z k z_k zk，再通过 z k z_k zk 推导得到原先的 θ k \theta_k θk。
如何从 2 M − 2 2M-2 2M−2 个根中确定 K K K 个根
那么如何从 2 M − 2 2M-2 2M−2 个根中确定 K K K 个根？这个问题大部分的论文和博客都一笔带过了。从前面的讨论可知，多项式系数是由 G \mathbf{G} G 的多条对角线求和得到，同时 G \mathbf{G} G 是 Hermitian 矩阵，因此以下式子可以得到：
a p = a − p ∗ a_p = a_{-p}^* ap=a−p∗
这个等式意味着在 2 M − 1 2M-1 2M−1 个系数 { a p : p = − M + 1 , ⋯ , M − 1 } \{a_p:p=-M+1,\cdots,M-1\} {ap:p=−M+1,⋯,M−1} 中，前 M − 1 M-1 M−1 个和后 M − 1 M-1 M−1 个系数是前后共轭对称，同时正中间的系数是实数。
我们继续假设 M = 3 M=3 M=3， P − 1 ( z ) = 0 P^{-1}(z)=0 P−1(z)=0 可以进一步表示如下：
P − 1 ( z ) = a 2 z − 2 + a 1 z − 1 + a 0 + a 1 ∗ z 1 + a 2 ∗ z 2 = 0 P^{-1}(z) = a_{2}z^{-2}+a_{1}z^{-1} + a_{0} + a_{1}^*z^{1} + a_{2}^*z^{2}=0 P−1(z)=a2z−2+a1z−1+a0+a1∗z1+a2∗z2=0
如此我们分析 P − 1 ( 1 / z ∗ ) P^{-1}(1/z^*) P−1(1/z∗)，可得：
P − 1 ( 1 / z ∗ ) = a 2 ( z ∗ ) 2 + a 1 ( z ∗ ) 1 + a 0 + a 1 ∗ ( z ∗ ) − 1 + a 2 ∗ ( z ∗ ) − 2 = [ P − 1 ( z ) ] ∗ = P − 1 ( z ) = 0 \begin{aligned} P^{-1}(1/z^*) &= a_{2}(z^*)^{2}+a_{1}(z^*)^{1} + a_{0} + a_{1}^*(z^*)^{-1} + a_{2}^*(z^*)^{-2} \\ &=[P^{-1}(z)]^* = P^{-1}(z) = 0 \end{aligned} P−1(1/z∗)=a2(z∗)2+a1(z∗)1+a0+a1∗(z∗)−1+a2∗(z∗)−2=[P−1(z)]∗=P−1(z)=0
这意味着假若 z 1 = ρ e j φ z_1 = \rho e^{\mathrm{j}\varphi} z1=ρejφ 是 P − 1 ( z ) = 0 P^{-1}(z)=0 P−1(z)=0 的根，那么 z 2 = 1 / z 1 ∗ = 1 / ρ e j φ z_2 = 1/z_1^* = 1/\rho e^{\mathrm{j}\varphi} z2=1/z1∗=1/ρejφ 同样是 P − 1 ( z ) = 0 P^{-1}(z)=0 P−1(z)=0 的根。观察 z 1 z_1 z1 和 z 2 z_2 z2 在复平面的位置，将会观察得到 z 1 z_1 z1 和 z 2 z_2 z2 是关于单位圆有一个类似对称的关系；简单来说，这个现象是因为 z 1 z_1 z1 和 z 2 z_2 z2 是幅值互为倒数而相位相等的关系，因此它们就像是关于 e j φ e^{\mathrm{j}\varphi} ejφ 对称一样（ e j φ e^{\mathrm{j}\varphi} ejφ 的幅值是 1 1 1，因此该点在单位圆上）。
综上所述，通过 P − 1 ( z ) P^{-1}(z) P−1(z) 得到 2 M − 2 2M-2 2M−2 个根，它们是关于单位圆对称的 M − 1 M-1 M−1 对根，因此一定有 K K K 对根在单位圆附近，所以我们只需要从 2 M − 2 2M-2 2M−2 个根中找 M − 1 M-1 M−1 个处于单位圆内的根（找 M − 1 M-1 M−1 个处于单位圆外的根同样是可以的，因为角度信息其实只存在于 z k z_k zk 的相位中，与幅值无关），最后确定最接近单位圆的 K K K 个根就可以确定 z k z_k zk。
从复数域上观察 2 M − 2 2M-2 2M−2 个根的分布
我们从实验中进一步观察 2 M − 2 2M-2 2M−2 个根的分布，matlab 代码实现如下：
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 %% Signals S = exp(1j*2*pi*randn(K,T)); % signal matrix A = exp(-1j*(0:M-1)'*2*pi*d/lambda*sind(theta)); % steering vector matrix N = noise_flag.*sqrt(10.^(-SNR/10))*(1/sqrt(2))*(randn(M,T)+1j*randn(M,T)); % noise matrix X = A*S+N; % received matrix R = X*X'/T; % covariance matrix %% DoA:root-MUSIC [U,~] = svd(R); % SVD Un = U(:, K+1:end); % noise subspace matrix Gn = Un*Un'; coe = arrayfun(@(i) sum(diag(Gn, M-i)),(1:2*M-1)); r = roots(coe); % 2M-2 roots %% plot dis = sort(abs(r)-1); disp(dis); cnt = sum(dis
