创建网站的英语,石家庄微信网站建设,wordpress显示时间代码,建筑学是十大最烂专业之一吗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 算法#xff0c;需要理解多项式求… 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)a0a1x⋯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:i0,⋯,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:k1,⋯,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} zejω我们将会得到 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]Ta(θ)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^{-M1}] \mathbf{G} [1,z,z^{2},\cdots,z^{M-1}]^T \\ \sum_{m 0}^{M-1} \sum_{n0}^{M-1} z^{-m} \mathbf{G}_{[m,n]} z^{n} \\ \sum_{m 0}^{M-1} \sum_{n0}^{M-1} z^{n-m} \mathbf{G}_{[m,n]} \\ \sum_{p-M1}^{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−M1]G[1,z,z2,⋯,zM−1]Tm0∑M−1n0∑M−1z−mG[m,n]znm0∑M−1n0∑M−1zn−mG[m,n]p−M1∑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−np∑G[m,n] 到这里我们已经可以看出传统 MUSIC 算法对 P ( θ ) P(\theta) P(θ) 求峰值其实等价于对 P − 1 ( z ) P^{-1}(z) P−1(z) 求根为了方便大家的理解我们令 M 3 M3 M3此时会得到一条简单的式子 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−2a1z−1a0z0a−1z1a−2z2 可以看出其实 P − 1 ( θ ) P^{-1}(\theta) P−1(θ) 是一个 2 M − 1 5 2M-1 5 2M−15 项的多项式但还存在一个问题 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−10 将会得到 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−10 求根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}^* apa−p∗ 这个等式意味着在 2 M − 1 2M-1 2M−1 个系数 { a p : p − M 1 , ⋯ , M − 1 } \{a_p:p-M1,\cdots,M-1\} {ap:p−M1,⋯,M−1} 中前 M − 1 M-1 M−1 个和后 M − 1 M-1 M−1 个系数是前后共轭对称同时正中间的系数是实数。 我们继续假设 M 3 M3 M3 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−2a1z−1a0a1∗z1a2∗z20 如此我们分析 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∗)2a1(z∗)1a0a1∗(z∗)−1a2∗(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} z21/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*SN; % received matrix
R X*X/T; % covariance matrix%% DoA:root-MUSIC
[U,~] svd(R); % SVD
Un U(:, K1: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(dis0);
disp(cnt); % 记录单位圆内的根个数% 提取实部和虚部
realPart real(r);
imaginaryPart imag(r);% 绘制复平面
figure;
scatter(realPart, imaginaryPart, filled);
hold on;% 绘制单位圆
theta linspace(0, 2*pi, 100);
unitCircleReal cos(theta);
unitCircleImag sin(theta);
plot(unitCircleReal, unitCircleImag, r--, LineWidth, 1);xlabel(实部);
ylabel(虚部);
title(复平面上的复数点和单位圆);
grid on;
box on;%% find K roots
r r(abs(r)1);
[~, idx] sort(abs(abs(r)-1));
z angle(r(idx));
theta sort(asin(-z(1:K)/2/pi/d*lambda)/pi*180).;设 M 16 M16 M16 和 K 2 K2 K2根的分布如下图所示可以看到 2 M − 2 30 2M-2 30 2M−230 个根其中 2 K 4 2K4 2K4 个接近单位圆的根对应着估计角度
