root MUSIC 算法补充说明


  这篇笔记是上一篇关于 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\} {aiC: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,ej2πdsinθ/λ,,ej2π(M1)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 P1(θ)=0
  在接下来的讨论中,我们令 P − 1 ( θ ) = a H ( θ ) G a ( θ ) P^{-1}(\theta) = \mathbf{a}^H(\theta)\mathbf{G}\mathbf{a}(\theta) P1(θ)=aH(θ)Ga(θ),此时我们可以知道,MUSIC 算法满足 G ≜ U n U n H \mathbf{G} \triangleq \mathbf{U}_n\mathbf{U}^H_n GUnUnH,而 Capon 算法满足 G ≜ R − 1 \mathbf{G} \triangleq \mathbf{R}^{-1} GR1。需要注意的是无论是 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)P1(z)=[1,z,z2,,zM1]T=a(θ)=aH(z)Ga(z)=P1(θ)
接下来我们展开 P − 1 ( z ) P^{-1}(z) P1(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} P1(z)=aH(z)Ga(z)=[1,z,(z)2,,(z)M1]G[1,z,z2,,zM1]T=[1,z1,z2,,zM+1]G[1,z,z2,,zM1]T=m=0M1n=0M1zmG[m,n]zn=m=0M1n=0M1znmG[m,n]=p=M+1M1apzp
其中 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]} apmn=pG[m,n]
  到这里我们已经可以看出,传统 MUSIC 算法对 P ( θ ) P(\theta) P(θ) 求峰值,其实等价于对 P − 1 ( z ) P^{-1}(z) P1(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} P1(z)=a2z2+a1z1+a0z0+a1z1+a2z2
可以看出,其实 P − 1 ( θ ) P^{-1}(\theta) P1(θ) 是一个 2 M − 1 = 5 2M-1 = 5 2M1=5 项的多项式,但还存在一个问题, P − 1 ( θ ) P^{-1}(\theta) P1(θ) 中存在负整数次数,我们令 P − 1 ( z ) z M − 1 P^{-1}(z)z^{M-1} P1(z)zM1 将负整数次数消除即可,操作前后,求根的结果是一样的,因此我们可以说 P − 1 ( z ) z M − 1 P^{-1}(z)z^{M-1} P1(z)zM1 是一个一元的 2 M − 1 2M-1 2M1 项的 2 M − 2 2M-2 2M2 次的多项式。更进一步地,我们可以说求解 P − 1 ( z ) z M − 1 = 0 P^{-1}(z)z^{M-1}=0 P1(z)zM1=0 将会得到 2 M − 2 2M-2 2M2 个根,从已知条件我们知道,其中 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 P1(z)zM1=0 求根,root MUSIC 算法所做的,就是利用 G \mathbf{G} G 的多条对角线求和得到对应的多项式系数,从而求解得 2 M − 2 2M-2 2M2 个根,接着筛选得到合适的 K K K 个根 z k z_k zk,再通过 z k z_k zk 推导得到原先的 θ k \theta_k θk

如何从 2 M − 2 2M-2 2M2 个根中确定 K K K 个根

  那么如何从 2 M − 2 2M-2 2M2 个根中确定 K K K 个根?这个问题大部分的论文和博客都一笔带过了。从前面的讨论可知,多项式系数是由 G \mathbf{G} G 的多条对角线求和得到,同时 G \mathbf{G} G 是 Hermitian 矩阵,因此以下式子可以得到:
a p = a − p ∗ a_p = a_{-p}^* ap=ap
这个等式意味着在 2 M − 1 2M-1 2M1 个系数 { a p : p = − M + 1 , ⋯ , M − 1 } \{a_p:p=-M+1,\cdots,M-1\} {ap:p=M+1,,M1} 中,前 M − 1 M-1 M1 个和后 M − 1 M-1 M1 个系数是前后共轭对称,同时正中间的系数是实数。
  我们继续假设 M = 3 M=3 M=3 P − 1 ( z ) = 0 P^{-1}(z)=0 P1(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 P1(z)=a2z2+a1z1+a0+a1z1+a2z2=0
  如此我们分析 P − 1 ( 1 / z ∗ ) P^{-1}(1/z^*) P1(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} P1(1/z)=a2(z)2+a1(z)1+a0+a1(z)1+a2(z)2=[P1(z)]=P1(z)=0
这意味着假若 z 1 = ρ e j φ z_1 = \rho e^{\mathrm{j}\varphi} z1=ρejφ P − 1 ( z ) = 0 P^{-1}(z)=0 P1(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 P1(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) P1(z) 得到 2 M − 2 2M-2 2M2 个根,它们是关于单位圆对称的 M − 1 M-1 M1 对根,因此一定有 K K K 对根在单位圆附近,所以我们只需要从 2 M − 2 2M-2 2M2 个根中找 M − 1 M-1 M1 个处于单位圆内的根(找 M − 1 M-1 M1 个处于单位圆外的根同样是可以的,因为角度信息其实只存在于 z k z_k zk 的相位中,与幅值无关),最后确定最接近单位圆的 K K K 个根就可以确定 z k z_k zk

从复数域上观察 2 M − 2 2M-2 2M2 个根的分布

  我们从实验中进一步观察 2 M − 2 2M-2 2M2 个根的分布,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<0);
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 M=16 M=16 K = 2 K=2 K=2,根的分布如下图所示,可以看到 2 M − 2 = 30 2M-2 = 30 2M2=30 个根,其中 2 K = 4 2K=4 2K=4 个接近单位圆的根对应着估计角度:
复平面上的复数点和单位圆

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.hqwc.cn/news/474480.html

如若内容造成侵权/违法违规/事实不符,请联系编程知识网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

浅析太阳能电池量子效率测试系统的主要组成部分

太阳能电池量子效率测试系统是用于对太阳能电池进行量子效率测试的设备。量子效率是指太阳能电池在接收光照射时&#xff0c;将光子转化为电子的效率。太阳能电池的量子效率越高&#xff0c;其转化光能为电能的效率就越高。主要由以下几个组成部分构成&#xff1a; 光源&#x…

阿里云配置服务器详细指南_2024年CPU内存带宽配置选择

阿里云服务器配置怎么选择&#xff1f;根据实际使用场景选择&#xff0c;个人搭建网站可选2核2G配置&#xff0c;访问量大的话可以选择2核4G配置&#xff0c;企业部署Java、Python等开发环境可以选择2核8G配置&#xff0c;企业数据库、Web应用或APP可以选择4核8G配置或4核16G配…

【PyQt】在PyQt5的界面上集成matplotlib绘制的图像

文章目录 0 前期教程1 概述2 matplotlib2.1 库导入2.2 图片的各个部分解释2.3 代码风格2.4 后端 3 集成matplotlib图像到pyqt界面中3.1 使用到的模块3.2 理解Qt Designer中的“控件提升”3.3 界面与逻辑分离的思路3.4 扩展 0 前期教程 【PyQt】PyQt5进阶——串口上位机及实时数…

用例是独立的-《软件方法》自测题解析38

DDD领域驱动设计批评文集 做强化自测题获得“软件方法建模师”称号 《软件方法》各章合集 《软件方法》第5章自测题2 8 [ 多选题 ] 根据以下业务序列图&#xff0c;请问属于“一卡通系统”用例的有&#xff1a;  A) 外来办事人员→登记  B) 一卡通系统→判断黑名单 …

GIS工具maptalks——初始化

GIS工具maptalks开发手册(一)——hello world初始化 为何使用maptalks&#xff1f; ​ Maptalks项目是一个HTML5的地图引擎, 基于原生ES6、Javascript开发的二三维一体化地图。 通过二维地图的旋转、倾斜增加三维视角&#xff0c;通过插件化设计, 能与其他图形库echarts、d3.j…

如何在 Linux 中安装 s3cmd 并管理 Amazon s3 存储桶

S3&#xff0c; – 简单存储服务- 是亚马逊的存储服务&#xff0c;为 IT 团队提供一种安全、可扩展且可靠的方式来存储和检索云上的文件和文件夹。 S3 可确保数据在需要时可用并随着需求的增长而扩展&#xff0c;从而帮助您充分利用数据。 通常&#xff0c;在登录到您的 AWS 账…

Mybatis源码分析

Mybatis源码分析 第一章、回顾 1. 课程中工具的版本 1. JDK8 2. IDEA2018.3 3. Maven3.5.3 4. MySQL 5.1.48 --> MySQL 5Mybatis 3.4.62. Mybatis开发的简单回顾 1. Mybatis做什么&#xff1f;Mybatis是一个ORM类型框架&#xff0c;解决的数据库访问和操作的问题&#xf…

Manifest merger failed with multiple errors, see logs

问题 Manifest merger failed with multiple errors, see logs详细问题 笔者进行Android 项目开发&#xff0c;修改AndroidManifest.xml代码后&#xff0c;控制台报错 AndroidManifest.xml报错核心代码 <manifest><uses-permission android:name"android.perm…

[03] Vue指令(2)

目录 属性绑定指令(v-bind)列表渲染指令(v-for)v-for中的key 双向绑定指令(v-model) 属性绑定指令(v-bind) 作用&#xff1a;动态设置html的标签属性&#xff0c;比如&#xff1a;src、url、title语法&#xff1a;v-bind:属性名“表达式”v-bind:可以简写成 > : 比如&…

C语言函数调用的流程

函数调用时&#xff0c;进程的上下文会切换到被调函数&#xff0c;当被调函数执行完毕之后再切换回去 函数调用时代码的执行流程

电子商务跨境电商大数据的关键技术之—主流电商大数据采集

大数据采集是指通过各种技术手段和工具收集、获取和提取大规模数据的过程。在信息时代&#xff0c;各种互联网、物联网、移动设备等的普及和应用&#xff0c;产生了海量的数据&#xff0c;这些数据被称为大数据。大数据采集就是对这些数据进行收集和抓取&#xff0c;以获得有意…

SpringBoot+Vue项目部署上线

部署前准备 注册京东云 京东云: https://www.jdcloud.com/ 117.72.32.65vue本地部署 新建文件.env.development VUE_APP_BASEURLhttp://localhost:9191新建文件 .env.production VUE_APP_BASEURLhttp://117.72.32.65:9191main.js 设置全局变量$baseUrl Vue.prototype.$baseUrl…