阵列信号处理-波达方向DOA-子空间方法

• 前言
• 子空间和子空间数据模型
• • 子空间
  • 子空间数据模型
  • 子空间模型的参数估计
  • • 1.参数建模
    • 信号采样
    • 参数估计
  • DOA 估计问题
  • • Deterministic maximum likelihood方法
    • Beam forming方法
• MUSIC和ESPRIT
• • MUSIC
• SRP-PHAT
• CSSM
• WAVES
• TOPS
• FRIDA
• 术语和约定-notational conventions
• • ARMA
  • FSD
  • low-rank && full-rank
  • azimuth and elevation angles
  • wavefronts：planar or curvature
  • complex：in-phase and quadrature
  • coherent & incoherent
  • covariance协方差
  • 协方差矩阵
  • 特征向量、特征值
  • Tr. 矩阵的迹
  • Maximum likelihood estimation
  • Singular Value Decomposition（奇异值分解）
• 参考文献

前言

波达方向-Direction Of Arrival是研究波束形成的重要课题，引用之前的老图，DOA要估算出来的就是两个角度：俯仰角 φ \varphi φ和方位角 θ \theta θ，而如果是全向麦克风组成的线性阵列（z轴），那么方位角就可以省略了，只研究俯仰角就可以了，所以很多算法简化假设条件，然而实际中无法省略，不过可以通过阵列的摆放，设计的估算角度范围小一些，算法也会容易一些。

经典的DOA估算方法有波束形成测向方法（base），Capon最小功率估计器，ML极大似然估计器，MUSIC 多重信号分类方法，ESPRIT旋转不变量信号参数估计方法，随着时代进步，开源代码提供了丰富晚上的功能，本文跟随pyroomacoustic实现的方法，在参考文献中找出来这些神作，学习记录，以备后用。要研究这些算法，我们必须先了解一下子空间的基础知识。

子空间和子空间数据模型

说句题外话，我已经记不起来大学里学没学过这个概念了，不过在科幻电影主导荧幕的21世纪，提到多维空间的子空间是神秘而且令人兴奋的。只是线代里的子空间定义非常滴扼要，以至于看完了概念仍然一头雾水，尽管如此还是现学现卖了抄写一下公式吧。

子空间

【from 百度百科】假设设W是域P上的线性空间V的一个非空子集合，若对于V中的加法及域P与V的纯量乘法构成域P上的一个线性空间，则称W为V的线性子空间（或向量子空间），或简称子空间。
换句话讲，V的非空子集W是子空间的充分必要条件是：
（1）子集合W的任意两个向量α与β之和α+β仍是W中的向量；
（2）域P的任一数k与子集合W的任意一个向量α的积kα仍是W中的向量。
用直观的概念来理解，就是三维坐标中，任意由两个轴组成的平面即是这三维空间的二维子空间，所以多维空间中，子空间也有很多性质，都可以用三维空间来联想一下，方便记忆和理解。

子空间数据模型

【9】将阵列信号看成矩阵 A A A的参数估计问题，假设：
x ( t ) = A ( η ) s ( t ) + n ( t ) x(t)= \mathbf{A}(\boldsymbol{\eta})s(t)+n(t) x(t)=A(η)s(t)+n(t)
这里先重新声明一下变量的意义和所属集合：
x ( t ) ∈ C M x(t)\in \Bbb{C}^M x(t)∈CM，可以理解为有M个传感器变换的矩阵，在 t t t时刻采集到的信号；
s ( t ) ∈ C k s(t)\in \Bbb{C}^k s(t)∈Ck，是远场输入的k个信号， n ( t ) ∈ C M n(t)\in \Bbb{C}^M n(t)∈CM是加性噪声；
A ∈ C M × k \mathbf{A}\in \Bbb{C}^{M\times k} A∈CM×k是转换矩阵，假设在观测窗内具有时不变特性，而 η ∈ R q \boldsymbol{\eta}\in \Bbb{R}^q η∈Rq是决定这个矩阵的参数向量，如果是线阵理想的情况，就是入射角 θ \theta θ。
另外分析此类问题，还是离不开假设：

1. A ( η ) \mathbf{A}(\boldsymbol{\eta}) A(η)在观测周期时不变；
2. 模型在 A \mathbf{A} A和 s s s上满足双线性（bilinear）；
3. 加性噪声；

基于以上所有的假设和定义，那么当 k < M k<M k<M，或者说入射信号数量少于观测传感器数量的情况，那么信号 x ( t ) ∈ C M x(t)\in \Bbb{C}^M x(t)∈CM能够被压缩到 C M \Bbb{C}^M CM一个 k k k维子空间 C k \Bbb{C}^k Ck，这个子空间被称为“信号子空间”。而噪声被假设为在所有维度上都有能量分布的，所以，此时阵列信号公式可以理解为满秩（full-rank）噪声数据模型所污染的低秩（low-rank）信号，这就是基于子空间的数据模型定义。
这个模型的几何意义，可以用三维坐标来说明，下图引自知乎上的博文：

如果 x ( t ) ∈ C 3 x(t)\in \Bbb{C}^3 x(t)∈C3被变换到子空间 C 2 \Bbb{C}^2 C2子空间来分析，那么弥散在整个三维空间的噪声大部分被降维规避了，分布于 C 2 \Bbb{C}^2 C2上的噪声能量势必远小于整个三维空间。

子空间模型的参数估计

本段基本翻译理解【9】，基于子空间模型的参数估计大致分三步：

1.参数建模

\quad 建模一个参数 η \boldsymbol{\eta} η的估计向量 η = [ η 1 , . . . , η k ] \boldsymbol{\eta}=[\boldsymbol{\eta}_1,...,\boldsymbol{\eta}_k] η=[η1​,...,ηk​]由（列）向量组成（矩阵） A ( η ) = [ a ( η 1 ) , . . . , a ( η k ) ] \mathbf{A}(\boldsymbol{\eta})=[\boldsymbol{a(\eta_1}),..., \boldsymbol{a(\eta_k})] A(η)=[a(η1​),...,a(ηk​)] A ( η ) \mathbf{A}(\boldsymbol{\eta}) A(η)中的列向量 a ( η k ) \boldsymbol{a(\eta_k}) a(ηk​)表示为第 k k k个信号源的响应，线阵的方向列向量 a ( θ k ​ ) = [ 1 , e − j 2 π d s i n θ k ​ / λ 0 ​ , . . e − j 2 π ( M − 1 ) d s i n θ k ​ / λ 0 ​ ] a(\theta_k​)=[1,e^{−j2πdsinθk​/λ0}​,..e^{−j2π(M−1)dsinθk​/λ0​}] a(θk​​)=[1,e−j2πdsinθk​/λ0​,..e−j2π(M−1)dsinθk​/λ0​]就是在线阵DOA估计中的一个实例。而一般情况下， η k \boldsymbol{\eta_k} ηk​本身也是多维向量，有可能包含坐标信息，信号载波频率等等，从算法上只要这个向量的维度（假设 p p p）小于阵元数量 M M M，那么 a ( η k ) \boldsymbol{a(\eta_k}) a(ηk​)能够在 C M \Bbb{C}^M CM空间中，以 η k \boldsymbol{\eta_k} ηk​作为变量，绘制一个 p p p维曲面，这个曲面被叫做阵列流形。数学上定义如下为 A = { a ( η k ) : η ∈ Θ } \mathscr{A}=\{\boldsymbol{a(\eta_k}): \boldsymbol{\eta} \in \boldsymbol{\Theta}\} A={a(ηk​):η∈Θ}这里 Θ \boldsymbol{\Theta} Θ是参数向量所有可能取值的集合。很显然线阵的DOA情况 p = 1 p=1 p=1， η k = θ k \boldsymbol{\eta_k}=\theta_k ηk​=θk​，那么 A \mathscr{A} A是一个一维的曲线，这个曲线向绳子一样贯穿整个 C M \Bbb{C}^M CM空间，在【9】中有一个图示如下：

从几何上来看，窄带DOA估计问题可以理解为流形和信号子空间的交汇的地方（懵懂中），这就是我们期望的估计点。

一般情况下，如果 A \mathscr{A} A是确定的（unambiguous），并且阵元的个数 M M M多于信源的个数 k k k,那么矩阵 A ( η ) \mathbf{A}(\boldsymbol{\eta}) A(η)将是一个k维满秩的矩阵。集合阵元 X \mathit{X} X和信源 S \mathit{S} S的信号向量可以定义： X = A ( η ) S \mathit{X}=\mathbf{A}(\boldsymbol{\eta})\mathit{S} X=A(η)S按照子空间的理论，这个瞬间观察值（snapshot） 是矩阵 A ( η ) \mathbf{A}(\boldsymbol{\eta}) A(η)列的线性组合，或者说是每次观测都被限定在 C M \Bbb{C}^M CM的 k k k维子空间，这个子空间是被 A ( η ) \mathbf{A}(\boldsymbol{\eta}) A(η)的 k k k列所定义的。而如果信源 S \mathit{S} S本身是满秩的（ k k k），那么 A ( η ) \mathbf{A}(\boldsymbol{\eta}) A(η)的生成子空间将等价于 A \mathbf{A} A的生成子空间 s p a n ( A ) = s p a n ( A ( η ) ) span(\mathbf{A})=span(\mathbf{A}(\boldsymbol{\eta})) span(A)=span(A(η))，那么对阵元的观察向量 X \mathit{X} X将填满这个低秩子空间，为后面的矩阵参数估计提供了数学假设和基础。在加入了零均值加性噪声之后的阵元信号 X = A ( η ) S + N \mathit{X}=\mathbf{A}(\boldsymbol{\eta})\mathit{S}+\mathit{N} X=A(η)S+N

噪声就像鬼魂一样无处不在，看不见摸不着。但如鬼最好画，噪声也被各路大神描绘的惟妙惟肖，在这里将噪声假设成复稳态循环高斯随机过程，采样点之间完全不相关，以及空间协方差矩阵被定义为 σ 2 Σ = ξ { n ( t ) n ∗ ( t ) } \sigma^2\Sigma=\boldsymbol{\xi}\{n(t)n^*(t)\} σ2Σ=ξ{n(t)n∗(t)} ，这样来表示噪声的数学期望， Σ \Sigma Σ是归一化的矩阵即行列式 d e t ( Σ ) = 1 det(\Sigma)=1 det(Σ)=1， σ 2 \sigma^2 σ2表示的噪声能量。更大胆的假设是空间白噪声，那么用单位阵 I I I直接替代 Σ \Sigma Σ。
\quad 开完了噪声的玩笑，但引入噪声后的实际信号向量会变得满秩，就是无法用低秩序列来直接表征了。如果说不考虑噪声的情况，利用 N ≥ k N\geq k N≥k个数据向量所定义的信号子空间，寻找阵列流形和子空间的焦点作为解决方案。当噪声进来以后，就是要从数据中估计出来信号子空间，进而得到 η \eta η，让生成的阵列流形最佳适配这个估计。实值上，带噪信号的处理效率是子空间方法的关键问题。接下来我们看看带噪的子空间估计（SVD方法此处被提到），阵列协方差矩阵 R X X \boldsymbol{R}_{XX} RXX​可以被写成 R X X = E [ X X ∗ ] = A ( η ) R S S A ∗ ( η ) + σ 2 I \boldsymbol{R}_{XX}=\boldsymbol{E}[{XX^*}]=\mathbf{A}(\boldsymbol{\eta})\boldsymbol{R}_{SS}\mathbf{A}^*(\boldsymbol{\eta})+\sigma^2\mathbf{I} RXX​=E[XX∗]=A(η)RSS​A∗(η)+σ2I这里 R S S = E [ S S ∗ ] \boldsymbol{R}_{SS}=\boldsymbol{E}[{SS^*}] RSS​=E[SS∗]是信源协方差矩阵，这个矩阵是无法测量的，但从前面的推导，可以看出 A ( η ) R S S A ∗ ( η ) \mathbf{A}(\boldsymbol{\eta})\boldsymbol{R}_{SS}\mathbf{A}^*(\boldsymbol{\eta}) A(η)RSS​A∗(η)这个矩阵有一个低秩的特点，因为信源的个数也是茫然的，所以假设这是一个秩为 d ′ < M d^{\prime}<M d′<M的矩阵，通过特征值分解（eigendecomposition）得出 R X X = ∑ k = 1 M λ k e k e k ∗ \boldsymbol{R}_{XX}=\sum_{k=1}^M\lambda_ke_ke_k^* RXX​=k=1∑M​λk​ek​ek∗​
λ 1 > λ 2 > . . . . . λ M \quad\lambda_1>\lambda_2>.....\lambda_M λ1​>λ2​>.....λM​为特征值， { e k } \{e_k\} {ek​}是相对应的特征向量。根据 R X X = A ( η ) R S S A ∗ ( η ) + σ 2 I \boldsymbol{R}_{XX}=\mathbf{A}(\boldsymbol{\eta})\boldsymbol{R}_{SS}\mathbf{A}^*(\boldsymbol{\eta})+\sigma^2\mathbf{I} RXX​=A(η)RSS​A∗(η)+σ2I公式， R X X \boldsymbol{R}_{XX} RXX​是有一个秩为 d ′ d^{\prime} d′的矩阵再加上一个定标的单位阵，这样会让你容易的推算出最后的 M − d ′ M-d^\prime M−d′个特征值是非常小的，即 λ d ′ + 1 = λ d ′ + 2 = . . . . . λ M = σ 2 \lambda_{d^\prime+1}=\lambda_{d^\prime+2}=.....\lambda_M=\sigma^2 λd′+1​=λd′+2​=.....λM​=σ2，由此，定义 E S = [ e 1 , e 2 , . . . e d ′ ] \boldsymbol{E_S}=[e_1,e_2,...e_{d^{\prime}}] ES​=[e1​,e2​,...ed′​]为信号子空间，和定义 E N = [ e d ′ + 1 , e d ′ + 2 , . . . e M ] \boldsymbol{E_N}=[e_{d^{\prime}+1},e_{d^{\prime}+2},...e_M] EN​=[ed′+1​,ed′+2​,...eM​]作为噪声子空间，这两个空间是正教的（orthogonal complement）。而信号子空间包含于 A ( η ) \mathbf{A}(\boldsymbol{\eta}) A(η)的生成空间。到此，基于子空间估计的建模算是推导完毕。

信号采样

在测量空间，采样获取信号子空间的估计值 A ^ \hat\mathbf{A} A^；实际中我们是估计 R ^ X X \boldsymbol{\hat R}_{XX} R^XX​，那么基于N次采样的估计公式为 R ^ X X = 1 N ∑ t = 1 N x ( t ) x ∗ ( t ) \boldsymbol{\hat R}_{XX}=\frac{1}{N}\sum_{t=1}^Nx(t)x^*(t) R^XX​=N1​t=1∑N​x(t)x∗(t)经过特征分解后的 ( λ ^ k , e ^ k ) , 1 ≤ k ≤ M (\hat \lambda_k, \hat e_k), 1\leq k\leq M (λ^k​,e^k​),1≤k≤M是 R ^ X X \boldsymbol{\hat R}_{XX} R^XX​的特征对，其中对最大的特征值进行保存，形成了信号子空间的持续估计，而信号子空间的维度也是基于特征值分布得到的。经典的方法有似然比（ikelihood ratio），MDL(minimum description length)和AIC(akaike information criterion)。到此可以看出计算矩阵的特征值分解将是最大的负担。

参数估计

估计一组参数 η ^ \hat\boldsymbol{\eta} η^​，使得 A ( η ^ ) \mathbf{A}(\hat\boldsymbol{\eta}) A(η^​)能够从某个层面最好的匹配 A ^ \hat\mathbf{A} A^；但这个 η ^ \hat\boldsymbol{\eta} η^​如何与 A ^ \hat\mathbf{A} A^建立联系，最终得到一直想要却还没出现的DOA呢？

DOA 估计问题

了解了子空间的概念和一般的分析方法，这些方法可以用在估算信号的方向、波长、频率、幅度，甚至利用多普勒频移估算信号的速度。而用在DOA上，是如何建模的呢？这种方法叫做子空间拟合（subspace fitting）。【9】子空间拟合办法就是一个公式 [ A ^ , T ^ ] = arg ⁡ min ⁡ A , T ∣ ∣ M − A ( η ) T ∣ ∣ F 2 \left [\hat \mathbf A,\hat \mathbf T\right]=\arg \min\limits_{A,T} ||\mathbf M-\mathbf{A}(\boldsymbol{\eta}) \mathbf T||^2_F [A^,T^]=argA,Tmin​∣∣M−A(η)T∣∣F2​相比于之前的一堆公式，多出来的 M M M是一个 m × q m\times q m×q矩阵，这个矩阵要从采样数据中重建。 m × p m\times p m×p的矩阵是利用 η \eta η参数化的，联想一下线阵的角向量。 T \mathbf T T是一个任意 p × q p\times q p×q矩阵，其中一个办法是替换一个 T ^ = A † M \hat \mathbf T=\mathbf A^ \dagger \mathbf M T^=A†M，那个奇怪的符号表示conjugate transpose，又称埃尔米特共轭。用这个替换重新获得一个等价的问题 A ^ = arg ⁡ max ⁡ A T r { P A M M ∗ } , P A = A ( A ∗ A ) − 1 A ∗ \hat \mathbf A=\arg \max\limits_{A}Tr\{P_AMM^*\}, P_A=A(A^*A)^{-1}A^* A^=argAmax​Tr{PA​MM∗},PA​=A(A∗A)−1A∗ P A P_A PA​叫做投影矩阵，投影到 A \mathbf A A的列空间。
\quad 矩阵 M M M的选择和包含 A \mathbf A A的集合就成了解决子空间拟合的基本问题，在【9】中有一个概括的表格，比较了各种子空间拟合方法：
M & A { A ∈ A d } { A ∈ A } { A ∈ ξ } M M ∗ = R ^ X X D e t − M L B e a m f o r m i n g M L − E S P R I T M = E ^ S M D − M U S I C M U S I C T L S − E S P R I T M = E ^ S W o p t 1 / 2 W S F W e i g h t e d   M U S I C W e i g h t e d   E S P R I T \begin{array}{c|lcr} {M\&A} & \{A\in\mathscr{A}^d\} & \{A\in\mathscr{A}\} & \{A\in\boldsymbol{\xi}\}\\ \hline MM^*=\boldsymbol{\hat R}_{XX} & Det-ML & Beamforming & ML-ESPRIT \\ M=\hat E_S & MD-MUSIC & MUSIC & TLS-ESPRIT \\ M= \hat E_SW_{opt}^{1/2}& WSF & Weighted\ MUSIC & Weighted\ ESPRIT \end{array} M&AMM∗=R^XX​M=E^S​M=E^S​Wopt1/2​​{A∈Ad}Det−MLMD−MUSICWSF​{A∈A}BeamformingMUSICWeighted MUSIC​{A∈ξ}ML−ESPRITTLS−ESPRITWeighted ESPRIT​​

Deterministic maximum likelihood方法

这种方法不是基于统计模型设计的，这样数据阵 X N X_N XN​最大化似然函数被认为是这样一个等价问题 min ⁡ η , S N T r { ( X N − A ( η ) S N ) ∗ ( X N − A ( η ) S N ) } = min ⁡ η , S N ∣ ∣ ( X N − A ( η ) S N ) ∣ ∣ F 2 \min\limits_{\eta,S_N}Tr\{(X_N-\mathbf{A}(\boldsymbol{\eta}) \mathbf S_N)^*(X_N-\mathbf{A}(\boldsymbol{\eta}) \mathbf S_N)\}=\min\limits_{\eta,S_N}||(X_N-\mathbf{A}(\boldsymbol{\eta}) \mathbf S_N)||_F^2 η,SN​min​Tr{(XN​−A(η)SN​)∗(XN​−A(η)SN​)}=η,SN​min​∣∣(XN​−A(η)SN​)∣∣F2​一般认为 S N S_N SN​是确定的，只是未知，所以假设 S ^ N = A † X N \hat \mathbf S_N=\mathbf A^ \dagger \mathbf X_N S^N​=A†XN​。代入上式替换 S N \mathbf S_N SN​，经过推导（我寄几是推不出来）最后得出优化问题的公式 η ^ = arg ⁡ max ⁡ η T r { P A ( η ) R ^ X X } \hat \mathbf {\eta}=\arg \max\limits_{\eta}Tr\{ \boldsymbol{P_{A(\eta)}}\boldsymbol{\hat R}_{XX}\} η^​=argηmax​Tr{PA(η)​R^XX​}, R ^ X X = s a m p l e   c o v a r i a n c e   m a t r i x \boldsymbol{\hat R}_{XX}=sample\ covariance\ matrix R^XX​=sample covariance matrix。这个方法是上面表格的第一个方法，计算负担很重，实际用的已经不多。

Beam forming方法

这里是指传统的’delay-and-sum method’，核心思想就是寻找加权输出的最大功率角度，阵列会扫描所有可能的方向已获得最优解。从公式的角度是取建模归一化的输出功率 P ( η j ) P(\eta_j) P(ηj​)，利用列向量 a ( η j ) a(\eta_j) a(ηj​)以及矩阵方法求解 P ( η j ) = 1 N ∑ t = 1 N ∣ a ∗ ( η j ) x ( t ) ∣ 2 ∣ a ( η j ) ∣ 2 = a ∗ ( η j ) R ^ X X a ( η j ) a ∗ ( η j ) a ( η j ) = T r { P a ( η j ) R ^ X X } P(\eta_j)=\frac{1}{N}\sum_{t=1}^N\frac{|a^*(\eta_j)x(t)|^2}{|a(\eta_j)|^2}=\frac{a^*(\eta_j)\hat R_{XX}a(\eta_j)}{a^*(\eta_j)a(\eta_j)}=Tr\{ \boldsymbol{P_{a(\eta_j)}}\boldsymbol{\hat R}_{XX}\} P(ηj​)=N1​t=1∑N​∣a(ηj​)∣2∣a∗(ηj​)x(t)∣2​=a∗(ηj​)a(ηj​)a∗(ηj​)R^XX​a(ηj​)​=Tr{Pa(ηj​)​R^XX​}， P a ( η j ) = a ( a ∗ a ) − 1 a ∗ P_{a(\eta_j)}=a(a^*a)^{-1}a^* Pa(ηj​)​=a(a∗a)−1a∗，最大功率方向就是DOA的方向。比较DML方法，如果信源只有一个，两种方法是等价的；而分离度比较好的声源，两者的估计也相当。而要分辨靠得很近的声源辐射，阵列孔径的音素不得不考虑。

MUSIC和ESPRIT

Multiple Signal Classification多信号分类方法是最经典的空间谱估计方法，也是在加性噪声环境下成功的构建数据模型的第一个高精度算法。而Estimation of Signal Parameters via Rotational Invariance Techniques 则对MUSIC的精度和限制进行了改进，有人戏称ESPRIT是MUSIC的儿子【8】，两者都是以子空间（subspace）技术为基础的方法，认为信号和噪声是在不同的子空间。pyroomacoustic没有采用ESPRIT，猜想是因为MUSIC本身的稳定性和鲁棒性更好一些。

MUSIC

“音乐”方法是1981年Ralph Otto Schmidt在斯坦福的博士论文《A Signal Subspace Approach to Multiple Emmitter Location and Spectral Estimation》提出来的方法，距今已经有40年了。他在Multiple emitter location and signal parameter estimation文章里有算法的经典推导过程。

SRP-PHAT

Steered Response Power – Phase Transform

CSSM

Coherent Signal Subspace Method

WAVES

Weighted Average of Signal Subspaces

TOPS

Test of Orthogonality of Projected Subspaces

FRIDA

Finite Rate of Innovation Direction of Arrival

术语和约定-notational conventions

ARMA

Autoregressive moving average model 自回归滑动平均模型， 是研究时间序列的重要方法，由自回归模型（简称AR模型）与移动平均模型（简称MA模型）为基础“混合”构成。

FSD

fast subspace decomposition

low-rank && full-rank

azimuth and elevation angles

方位角和仰角

wavefronts：planar or curvature

complex：in-phase and quadrature

coherent & incoherent

covariance协方差

协方差（Covariance）在概率论和统计学中用于衡量两个变量的总体误差。而方差是协方差的一种特殊情况，即当两个变量是相同的情况。协方差表示的是两个变量的总体的误差，这与只表示一个变量误差的方差不同。 如果两个变量的变化趋势一致，也就是说如果其中一个大于自身的期望值，另外一个也大于自身的期望值，那么两个变量之间的协方差就是正值。 如果两个变量的变化趋势相反，即其中一个大于自身的期望值，另外一个却小于自身的期望值，那么两个变量之间的协方差就是负值。（摘自百度百科）

协方差矩阵

统计学与概率论中，协方差矩阵的每个元素是各个向量元素之间的协方差，是从标量随机变量到高维度随机向量的自然推广。

特征向量、特征值

Tr. 矩阵的迹

Maximum likelihood estimation

似然和概率是鸡和蛋的问题，在给定条件 θ \theta θ下计算事件出现的概率可以表达为 P ( x ∣ θ ) P(x|\theta) P(x∣θ)，而根据事件发生的概率（结果）反推出条件就是似然的过程 L ( θ ∣ x ) \mathcal{L}(\theta|x) L(θ∣x)，所以最大似然估计是利用已知的样本的结果，在使用某个模型的基础上，反推最有可能导致这样结果的模型参数值。很显然求解DOA的过程也可以利用最大似然估计的办法。

Singular Value Decomposition（奇异值分解）

参考文献

【1】.R. Schmidt, Multiple emitter location and signal parameter estimation, IEEE Trans. Antennas Propag., Vol. 34, Num. 3, pp 276–280, 1986
【2】.J. H. DiBiase, A high-accuracy, low-latency technique for talker localization in reverberant environments using microphone arrays, PHD Thesis, Brown University, 2000
【3】.H. Wang, M. Kaveh, Coherent signal-subspace processing for the detection and estimation of angles of arrival of multiple wide-band sources, IEEE Trans. Acoust., Speech, Signal Process., Vol. 33, Num. 4, pp 823–831, 1985
【4】.E. D. di Claudio, R. Parisi, WAVES: Weighted average of signal subspaces for robust wideband direction finding, IEEE Trans. Signal Process., Vol. 49, Num. 10, 2179–2191, 2001
【5】.Y. Yeo-Sun, L. M. Kaplan, J. H. McClellan, TOPS: New DOA estimator for wideband signals, IEEE Trans. Signal Process., Vol. 54, Num 6., pp 1977–1989, 2006
【6】.H. Pan, R. Scheibler, E. Bezzam, I. Dokmanic, and M. Vetterli, FRIDA: FRI-based DOA estimation for arbitrary array layouts, Proc. ICASSP, pp 3186-3190, 2017
【7】.ESPRIT-Estimation of Signal Parameters Via Rotational Invariance Techniques
【8】.ESPRIT Algorithm- The Son of MUSIC
【9】Subspace Methods for Directions-of-Arrival, A. Paulraj, B. Ottersten, R. Roy, A. Swindlehurst, G. Xu and T. Kailath
【10】Performance Analysis of MUSIC and ESPRIT DOA Estimation Algorithms for Adaptive Array Smart Antenna in MobileCommunication
【11】Direction-of-Arrival Estimation Methods: APerformance-Complexity Tradeoff Perspective
【】