这段时间做了一些和增强相关的任务，虽然我用的方法是NN-based方法，但是这套传统理论我看了之后还是很感兴趣。学了一些深蓝学院的课程之后，结合之前做的，大概梳理了一些东西出来。

这里可以假设声音属于平面波，方向为 a \textbf a a;假设有M个麦克风组成麦克风阵列 [ p 0 , . . . , p M − 1 ] [\textbf p_0,...,\textbf p_{M-1}] [p0​,...,pM−1​] 则接收到的信号为: x ( n , p ) = [ x ( n , p 0 ) x ( n , p 1 ) . . . x ( n , p M − 1 ) ] \textbf x(n,\textbf p) = \left[ \begin{array}{c} x(n,\textbf p_0) \\ x(n,\textbf p_1) \\ ... \\ x(n, \textbf p_{M-1}) \end{array} \right] x(n,p)=⎣⎢⎢⎡​x(n,p0​)x(n,p1​)...x(n,pM−1​)​⎦⎥⎥⎤​ 由于入射方向 a \textbf a a是已知的，就可以通过入射方向估计出信号到达每个麦克风的时延 τ \tau τ τ = a T p i c , i = 0 , . . . , N − 1 \tau = \frac{\textbf a^T\textbf p_i}{c}, i = 0,...,N-1 τ=caTpi​​,i=0,...,N−1 意为每个麦克风的坐标在入射方向上的投影距离，处以声速 c c c，如果规定 θ \theta θ是入射角与z轴的夹角， ϕ \phi ϕ为向量在x0y的投影与x的夹角，那么 a \textbf a a可以表示为 a = [ − s i n θ c o s ϕ − s i n θ s i n ϕ − c o s θ ] \textbf a = \left[\begin{array}{c} -sin\theta cos \phi \\ -sin\theta sin \phi \\ -cos\theta \end{array} \right] a=⎣⎡​−sinθcosϕ−sinθsinϕ−cosθ​⎦⎤​

由于实际接收到的信号有延时差异。那么，接收到的信号在时域上就可以表示为 x ( n , p ) = [ x ( n − τ 0 ) x ( n − τ 1 ) . . . x ( n − τ M − 1 ) ] \textbf x(n,\textbf p) = \left[ \begin{array}{c} x(n-\tau_0) \\ x(n-\tau_1) \\ ... \\ x(n-\tau_{M-1}) \end{array} \right] x(n,p)=⎣⎢⎢⎡​x(n−τ0​)x(n−τ1​)...x(n−τM−1​)​⎦⎥⎥⎤​ 时域上的延时可以直接在频域上的相位中体现出来 X ( ω ) = [ e − j ω τ 0 e − j ω τ 1 . . . e − j ω τ M − 1 ] X ( ω ) \textbf X(\omega) = \left[ \begin{array}{c} e^{-j\omega\tau_0}\\ e^{-j\omega\tau_1}\\ ...\\ e^{-j\omega\tau_{M-1}} \end{array} \right]X(\omega) X(ω)=⎣⎢⎢⎡​e−jωτ0​e−jωτ1​...e−jωτM−1​​⎦⎥⎥⎤​X(ω) ω τ i = ω a T p i c = ( w c a ) T p i = k T p i = 2 π λ a T p i \omega\tau_i = \omega \frac{\textbf a^T \textbf p_i}{c}= (\frac{w}{c}\textbf a)^T\textbf p_i = \textbf k^T \textbf p_i = \frac {2\pi}{\lambda}\textbf a^T\textbf p_i ωτi​=ωcaTpi​​=(cw​a)Tpi​=kTpi​=λ2π​aTpi​ 其中， λ \lambda λ是波长， k \textbf k k就被称为波数，是后面会经常用到的一个概念，因为它同时包含了声源方向和频率两个信息。则，上面的频域接收信号表达式可以表示为 X ( ω ) = [ e − j k T p 0 e − j k T p 1 . . . e − j k T p M − 1 ] X ( ω ) = V k ( k ) X ( ω ) \textbf X(\omega) = \left[ \begin{array}{c} e^{-j\textbf k^T\textbf p_0}\\ e^{-j\textbf k^T\textbf p_1}\\ ...\\ e^{-j\textbf k^T\textbf p_{M-1}} \end{array} \right]X(\omega) = \textbf V_k(\textbf k)X(\omega) X(ω)=⎣⎢⎢⎡​e−jkTp0​e−jkTp1​...e−jkTpM−1​​⎦⎥⎥⎤​X(ω)=Vk​(k)X(ω) 其中 V k ( k ) \textbf V_k(\textbf k) Vk​(k)被称为阵列流形矢量（array mainifold vector）

一个频率响应函数表示的是一个滤波器的频率相应，但是在阵列系统中，最关键的就是引入了空间的概念，也就是有频率-空间，两个维度。那么在求响应函数的时候也需要把空间的信息考虑进来。 Y ( ω ) = H T ( ω ) X ( ω ) = H T ( ω ) V k ( k ) X ( ω ) \textbf Y(\omega) = \textbf H^T(\omega)\textbf X(\omega) = \textbf H^T(\omega)\textbf V_k(\textbf k)X(\omega) Y(ω)=HT(ω)X(ω)=HT(ω)Vk​(k)X(ω) 其中 H ( ω ) \textbf H(\omega) H(ω)就是频率响应函数。那么就可以规定频率-波数响应（frequency-wavnumber response）函数为 γ ( ω , k ) = H T ( ω ) V k ( k ) \gamma(\omega,\textbf k) = \textbf H^T(\omega)\textbf V_k(\textbf k) γ(ω,k)=HT(ω)Vk​(k) 那我们经常听到另一个概念是beam pattern，意味波束方向图，实际上是频率-波数响应在球面 k = 2 π / λ k = 2\pi / \lambda k=2π/λ上进行取值 B ( ω : θ , ϕ ) = γ ( ω , k ) ∣ k = 2 π λ a B(\omega:\theta, \phi) = \gamma(\omega,\textbf k)|_{\textbf k = \frac{2\pi}{\lambda}\textbf a} B(ω:θ,ϕ)=γ(ω,k)∣k=λ2π​a​

ULA适用于麦克风阵列是一排的这种情况，这种就被称为线性阵列。假设所有的麦克风阵列均匀分布在z轴上，相邻两麦克风的间距为d，也就是每个麦克风的坐标就是 $p_{z_i} =(i - \frac {M-1}2)d, i=0,1,…,M-1 $ p x i = p y i = 0 p_{x_i} = p_{y_i} = 0 pxi​​=pyi​​=0 把坐标带入阵列流形矢量可以得到 V k ( k ) = [ e − j ( 0 − M − 1 2 ) k z d , e − j ( 1 − M − 1 2 ) k z d , . . . , e − j ( M − M − 1 2 ) k z d ] T \textbf V_k(\textbf k) =[e^{-j(0-\frac {M-1}2)k_zd},e^{-j(1-\frac {M-1}2)k_zd},...,e^{-j(M-\frac {M-1}2)k_zd}]^T Vk​(k)=[e−j(0−2M−1​)kz​d,e−j(1−2M−1​)kz​d,...,e−j(M−2M−1​)kz​d]T k z = − 2 π λ c o s θ k_z = -\frac {2\pi}\lambda cos\theta kz​=−λ2π​cosθ 则它的波束方向图是 B θ ( θ ) = H T V θ ( θ ) = e − j ( M − 1 2 ) 2 π d λ c o s θ ∑ i = 0 N − 1 H i e j i 2 π d λ c o s θ , 0 ≤ θ ≤ π B_\theta(\theta) = \textbf H^T\textbf V_\theta(\theta) = e^{-j(\frac {M-1}2)\frac{2\pi d}{\lambda}cos\theta}\sum_{i=0}^{N-1}H_i e^{ji\frac{2\pi d}\lambda cos\theta}, 0\le\theta\le\pi Bθ​(θ)=HTVθ​(θ)=e−j(2M−1​)λ2πd​cosθi=0∑N−1​Hi​ejiλ2πd​cosθ,0≤θ≤π 若 N = 11 , H i = 1 N , d = λ 2 N=11,H_i=\frac 1 N,d = \frac \lambda 2 N=11,Hi​=N1​,d=2λ​

可以看出，频率-波数响应函数和权重向量之间是傅立叶变化的关系 B ψ ( ψ ) = H T V ψ ( ψ ) = e − j ( M − 1 2 ) ψ ∑ i = 0 N − 1 H i e j i ψ = e − j ( M − 1 2 ) ψ ( ∑ i = 0 N − 1 H i e − j i ψ ) ∗ B_\psi(\psi) = \textbf H^T\textbf V_\psi(\psi) = e^{-j(\frac {M-1}2)\psi}\sum_{i=0}^{N-1}H_i e^{ji\psi}= e^{-j(\frac {M-1}2)\psi}(\sum_{i=0}^{N-1}H_i e^{-ji\psi})^* Bψ​(ψ)=HTVψ​(ψ)=e−j(2M−1​)ψi=0∑N−1​Hi​ejiψ=e−j(2M−1​)ψ(i=0∑N−1​Hi​e−jiψ)∗

如果我们希望 k = k T \textbf k = \textbf k_T k=kT​,那么频率-波数响应函数为 γ ( ω , k − k T ) \gamma(\omega,\textbf k -\textbf k_T) γ(ω,k−kT​)，阵列流形矢量 [ e − j ( k − k T ) T p 0 e − j ( k − k T ) T p 1 . . . e − j ( k − k T ) T p M − 1 ] \left[ \begin{array}{c} e^{-j(\textbf k-\textbf k_T)^T\textbf p_0}\\ e^{-j(\textbf k-\textbf k_T)^T\textbf p_1}\\ ...\\ e^{-j(\textbf k-\textbf k_T)^T\textbf p_{M-1}} \end{array} \right] ⎣⎢⎢⎡​e−j(k−kT​)Tp0​e−j(k−kT​)Tp1​...e−j(k−kT​)TpM−1​​⎦⎥⎥⎤​ 实际上也就可以视为，对权重信息 H H H(或者可以称为 w w w)，进行修改 w T = [ e − j k T p 0 e − j k T p 1 . . . e − j k T p M − 1 ] ⨀ w \textbf w_T = \left[ \begin{array}{c} e^{-j\textbf k^T\textbf p_0}\\ e^{-j\textbf k^T\textbf p_1}\\ ...\\ e^{-j\textbf k^T\textbf p_{M-1}} \end{array} \right]\bigodot \textbf w wT​=⎣⎢⎢⎡​e−jkTp0​e−jkTp1​...e−jkTpM−1​​⎦⎥⎥⎤​⨀w

主要是可以通过调整权重向量，使得旁瓣尽可能小 比如汉明窗 实际上我们所想要得到的一种情况，就是目标方向的响应很大，同时旁瓣很小，再有就是主瓣的宽度越小越好。这都是通过调整权重矩阵w来实现的。包括beamforming的一些方法。

假设知道了监听的方向，同时知道了干扰噪声的方向，那么我就会希望干扰噪声出现一个零点，也就是对干扰噪声方向不进行任何响应。着就用到的是零点调向技术。数学符号可以表达为 w T V k ( k T ) = 1 \textbf w^T\textbf V_k(\textbf k_T) =1 wTVk​(kT​)=1同时满足 w T V k ( k i ) = 0 , i = 1 , 2 , . . . , M 0 \textbf w^T\textbf V_k(\textbf k_i) =0, i=1,2,...,M_0 wTVk​(ki​)=0,i=1,2,...,M0​也就是有 M 0 M_0 M0​个方向我们希望没有响应。 一般来说这种条件数大于未知变量数的问题，我们只有一个最小二乘解，也就是对所有的条件都不一定满足，但是各个条件的误差平方和是最小的。然而这里我们明确希望的一点，一定要满足零点约束的话，问题就会随之改变。

于是这个问题可以近似为，已知一个期望权值向量 w d w_d wd​，这个向量一般来说我们可以就假设 w T w_T wT​为期望的权值向量，它满足了无畸变约束。现在要求一个向量 w w w，使得两者尽可能相近，同时满足零点约束的限制。于是公式可以写成 { m i n ∣ ∣ w d − w ∣ ∣ 2 s . t . w H C = 0 \left \{ \begin{aligned} &min ||\textbf w_d - \textbf w||^2 \\ &s.t. \textbf w^H\textbf C = \textbf 0 \end{aligned}\right . {​min∣∣wd​−w∣∣2s.t.wHC=0​ 其中 C \textbf C C是零点向量长成的空间。也就是零点列向量组成的矩阵。

这种问题的解就很经典，使用拉格朗日乘子法即可，这个方法后面也会经常用到。最后得到的结果是 w H = w d H ( I M − C ( C H C ) − 1 C H ) \textbf w^H = \textbf w_d^H(\textbf I_M - \textbf C(\textbf C^H \textbf C)^{-1}\textbf C^H) wH=wdH​(IM​−C(CHC)−1CH) 实际上可以很明显地看出来，最终的结果是 w d \textbf w_d wd​减去一个东西。I后面的一坨是一个投影矩阵，所以最终的结果就是 w d \textbf w_d wd​减去了 w d \textbf w_d wd​在零点约束空间 C \textbf C C上的投影。

投影矩阵： 满足“幂等矩阵”同时又是“共轭对称阵”，那么它就是投影矩阵 C ( C H C ) − 1 C H \textbf C(\textbf C^H\textbf C)^{-1}\textbf C^H C(CHC)−1CH和 I M − C ( C H C ) − 1 C H \textbf I_M - \textbf C(\textbf C^H\textbf C)^{-1}\textbf C^H IM​−C(CHC)−1CH都是投影矩阵，且两者相互正交

阵列麦克风信号处理实际上最关键的就是引入了空间这个概念，如何体现空间，就是通过信号到达麦克风的时间差体现的，时域的时间差可以直接等价地表现在频域的相移上。然后我们可以通过频域的窗函数去控制它的beam pattern，最终的目的是表现出一种期望方向的信号很大，其他方向或者干扰方向信号为0或者尽可能小的一种情况。

我们会经常听到这个名字，它的原理是最简单的，只要知道入射方向，就能得到信号到达每个麦克风的时间差，只要将这些麦克风的信号在时间上对齐即可。 所以它的做法就是，首先对每一个麦克风的信号进行延时操作。 y ′ ( n ) = 1 N ∑ i = 0 M − 1 x ( n + τ i ) y'(n) = \frac 1 N \sum_{i=0}^{M-1}x(n+\tau_i) y′(n)=N1​i=0∑M−1​x(n+τi​) 那么相应地，它的窗函数就是 h i ( n ) = 1 N δ ( n + τ i ) h_i(n) = \frac 1N \delta(n+\tau_i) hi​(n)=N1​δ(n+τi​)，频域表示为 H ( ω ) = 1 N e − j ω τ i H(\omega) = \frac 1Ne^{-j \omega\tau_i} H(ω)=N1​e−jωτi​

D-S 在很多地方都会使用，它虽然没有限制任何噪声方向，但是它很稳定。

这个方法的核心思想就是，找到权值向量 w \textbf w w，保持住目标方向响应不变，同时让让输出信号功率最小，公式表达为 { m i n w R x x w s . t . w H V k ( k s ) = 1 \left \{ \begin{aligned} &min \quad \textbf w \textbf R_{xx} \textbf w \\ &s.t. \textbf w^H\textbf V_k(\textbf k_s) = 1 \end{aligned}\right . {​minwRxx​ws.t.wHVk​(ks​)=1​ 这个方程的解法就是拉格朗日乘子法，得到的解为 w = R x x − 1 V k ( k s ) V k H ( k s ) R x x − 1 V k ( k s ) \textbf w = \frac {\textbf R_{xx}^{-1}\textbf V_k(\textbf k_s)}{\textbf V_k^H(\textbf k_s)\textbf R_{xx}^{-1}\textbf V_k(\textbf k_s)} w=VkH​(ks​)Rxx−1​Vk​(ks​)Rxx−1​Vk​(ks​)​ 保持目标方向不变，最小化输出功率就相当于最小化噪声功率，因此，也可以写成

w = R n n − 1 V k ( k s ) V k H ( k s ) R n n − 1 V k ( k s ) \textbf w = \frac {\textbf R_{nn}^{-1}\textbf V_k(\textbf k_s)}{\textbf V_k^H(\textbf k_s)\textbf R_{nn}^{-1}\textbf V_k(\textbf k_s)} w=VkH​(ks​)Rnn−1​Vk​(ks​)Rnn−1​Vk​(ks​)​ 这里面的 R n n \textbf R_{nn} Rnn​就是噪声自相关矩阵。 这个方法高度依赖于声源方向的估计，方向预测的准，效果就好。

我对GSC的理解就是，在DSB上减去了噪声项。给我的感觉有点像AEC的一些做法，不同的是，GSC的噪声是估计出来的。

可以直接来看一个GSC的实例，我觉得比理论还好理解一些。整个框架氛围三个模块：FBF、ABM、AIC。

其中会提及到语音存在概率这个东西，传统做法会有Minima Controlled Recursive Averaging (MCRA)方法。神经网络方法会有估计mask，对mask进行sigmoid去估计语音存在概率。

现在再从理论上来说说GSC，它的思想是把权重 w \textbf w w投影在两个相互正交的平面上，且其中一个平面是我们定义的（称作“约束子空间”），另一个平面是正交补空间，被称为“自适应子空间”。 w = w q − w p \textbf w = \textbf w_q - \textbf w_p w=wq​−wp​，其中 w q \textbf w_q wq​是约束子空间的投影， w p \textbf w_p wp​是自适应子空间的投影。

在上面的例子中， w q \textbf w_q wq​就是DSB的窗函数， w p w_p wp​就是ABM和AIC组合而成的。 w p = B ( B H B ) − 1 B H w = B w a w_p = \textbf B (\textbf B^H\textbf B)^{-1}\textbf B^H\textbf w = \textbf B \textbf w_a wp​=B(BHB)−1BHw=Bwa​具体来说，ABM对应了 B \textbf B B,AIC对应了 w a \textbf w_a wa​。两个变量都是通过自适应滤波器的方法求解得到的。

这个beamformer的目标函数为，最大化每个频点的信噪比 w G E V = a r g max ⁡ w ( k ) w H ( k ) R z z w ( k ) w H ( k ) R n n w ( k ) \textbf w_{GEV} = arg \max_{w(k)}\frac{\textbf w^H(k)\textbf R_{zz}\textbf w(k)}{\textbf w^H(k)\textbf R_{nn}\textbf w(k)} wGEV​=argw(k)max​wH(k)Rnn​w(k)wH(k)Rzz​w(k)​ 这里我们用 z z z表示干净语音信号。

这个问题实际上是一个广义瑞丽商问题。求解方法依旧是拉格朗日成子法 我们可以定义 R ( x ) = x T M x x T Q x R(x) = \frac {x^TMx}{x^TQx} R(x)=xTQxxTMx​，如果我们限制 x T Q x = 1 x^TQx = 1 xTQx=1，那么就可以用拉格朗日成子法来求极值。 结果是极点为 x = Q − 1 M x = Q^{-1}M x=Q−1M，对应的极值为 Q − 1 M Q^-1M Q−1M的特征值

最后解出来的结果为 w G E V ( ω ) = R n n − 1 ( ω ) R z z ( ω ) \textbf w_{GEV}(\omega)= \textbf R_{nn}^{-1}(\omega)\textbf R_{zz}(\omega) wGEV​(ω)=Rnn−1​(ω)Rzz​(ω)。

前面提到，阵列信号处理可以看成是空间、频率两个维度共同的响应。那么对于声源定位算法来说，就是查找适合的或者说是响应最大的空间方位。 所以很多空间定位算法，就是遍历一圈方向，找到最大的那个值，只不过是定义的值是不同的。

如果我们对DS的信号功率，在方向上进行遍历，选择一个最大的方向，这个算法就叫SRP算法。 k s = a r g max ⁡ k ∑ ω ∣ Y ( ω , k ) ∣ 2 \textbf k_s = arg \max_{k} \sum_\omega|Y(\omega,\textbf k)|^2 ks​=argkmax​ω∑​∣Y(ω,k)∣2

但是人们发现，更多时候，方向更多的时候是与相位相关的，与幅度无关，于是更倾向于使用SRP-PHANT的方法。由于两个麦克风的信号是高度相似的，所以可以对一个信号进行时间延迟（频域相移），计算相关函数，如果某个方向所有麦克风两两的相关响应最大，那么这个方向就是入射角。PHANT的方法使用的是广义互相关函数，也就是不考虑幅度大小。计算公式可以表示为： p ( k ) = ∑ i = 0 M − 1 ∑ l = 0 M − 1 ∑ ω X i ( ω ) X l ∗ ( ω ) ∣ X i ( ω ) X l ∗ ( ω ) ∣ e j ω τ i l p(\textbf k) = \sum_{i=0}^{M-1}\sum_{l=0}^{M-1}\sum_\omega\frac{X_i(\omega)X^*_l(\omega)}{|X_i(\omega)X^*_l(\omega)|}e^{j\omega\tau_{il}} p(k)=i=0∑M−1​l=0∑M−1​ω∑​∣Xi​(ω)Xl∗​(ω)∣Xi​(ω)Xl∗​(ω)​ejωτil​ 前面的 k \textbf k k体现在 τ i l \tau_{il} τil​上。 这个 τ i l \tau_{il} τil​就是TDOA。 τ i j = d c o s θ c \tau_{ij} = \frac{d cos\theta}{c} τij​=cdcosθ​ 还有一个有意思的理解，可以分享一下，公式前面是 X X ∗ ∣ X X ∗ ∣ \frac{XX^*}{|XX^*|} ∣XX∗∣XX∗​，实际上就是相位差（IPD），我们谈绝对相位是没有意义的，就像我们提绝对时间也是没有意义的是一样的，总是拿一个麦克风作为参考麦克风再去谈相对时延或者相位差，这才是有意义的。也就是前面是信号算出来的相位差，后面的 τ i l \tau_{il} τil​是我们估计的两个麦克风的到达时间差在相位的体现，去算没有归一化的cosine距离。

我们如果不用DS去算，而是使用MVDR去进行计算，那么得到的就是Capon谱 最优权重矩阵 w c \textbf w_c wc​下的阵列输出功率为 P ( k s ) = P c a p o n ( k ) = w c H R x x ( ω ) w c = 1 V k H ( k s ) R x x − 1 V k ( k s ) P(\textbf k_s) = P_{capon} (\textbf k)= \textbf w_c^H\textbf R_{xx}(\omega)\textbf w_c = \frac 1{\textbf V_{\textbf k}^H(\textbf k_s)\textbf R_{xx}^{-1}\textbf V_{\textbf k}(\textbf k_s)} P(ks​)=Pcapon​(k)=wcH​Rxx​(ω)wc​=VkH​(ks​)Rxx−1​Vk​(ks​)1​ k s = a r g max ⁡ k P c a p o n ( k ) \textbf k_s = arg \max_k P_{capon}(\textbf k) ks​=argkmax​Pcapon​(k)

现在在回过头来看看当时的google和amazon的做法，的确是很相似。google是让神经网络自己学一组steer vector，每个表示不同方向，在频域直接点乘。亚马逊是预设好几个方向的steer vector，这个是根据拓扑结构算出来的。后面就是方向选择问题，可以是拼一起，可以根据能量进行选择，或是attention。 在神经网络化的beamformer中，感觉大部分的神经网络的作用是用来预测speech或者是noise，或是预测T-F mask来得到speech和（或）noise。我们也看到了上面的很多公式中需要大量的speech或者是noise的自相关矩阵，这样估计出了mask之后就可以送给后面继续进行beamformer的估计。