简 介： 根据网络上一篇“基于广义互相关”的文章整理的所有广义互相关的公式。

关键词：广义互相关

在网络上有一篇关于基于广义互相关估计时间延时的PDF文件，是由Bhargava提交一篇综述报告，将用于估计时间延迟的广义互相关算法进行了总结。

■ 互相关估计时延

在空间上相距一定距离LLL的两个声音传感器获取远方声源发出的声音信号可以表示为：x1(t)=s1(t)+n1(t)x_1 \left( t \right) = s_1 \left( t \right) + n_1 \left( t \right)x1​(t)=s1​(t)+n1​(t)x2(t)=a⋅s1(t+D)+n2(t)x_2 \left( t \right) = a \cdot s_1 \left( {t + D} \right) + n_2 \left( t \right)x2​(t)=a⋅s1​(t+D)+n2​(t)

其中s1(t),n1(t),n2(t)s_1 \left( t \right),n_1 \left( t \right),n_2 \left( t \right)s1​(t),n1​(t),n2​(t)是实信号，可以看成联合平稳随机过程。信号s1(t)s_1 \left( t \right)s1​(t)与两个噪声n1(t),n2(t)n_1 \left( t \right),n_2 \left( t \right)n1​(t),n2​(t)之间是不相关的。

在很多应用中都基于估计时间延迟DDD。本文下面给出了一种最大似然（ML）估计算法来计算参数DDD，并与其他算法进行比较。

前面的模型是假设信号是联合平稳的，在实际应用中，通常信号的统计特性会发生缓慢变化。但只要在信号采集时间范围TTT之内保持近似平稳即可。

另外，估计参数时间延迟DDD也会随着时间变化发生缓慢变化。参数DDD的估计算子局限于在有限的观察时间范围TTT之内来估计DDD。其他的应用条件是关于信号和噪声的先验知识。通常，这些信号被忽略了。比如，在被动检测过程中，不像一些通讯应用中的问题，信号源的频谱是未知的，或者只是近似了解。

一种常用的计算时间延迟DDD的方法是通过计算两个传感器信号的互相关来获得的。基于时间延迟DDD，相对于传感器的信号波达方向也可以计算出来：Rx1,x2(τ)=E[x1(t)⋅x2(t)]R_{x_1 ,x_2 } \left( \tau \right) = E\left[ {x_1 \left( t \right) \cdot x_2 \left( t \right)} \right]Rx1​,x2​​(τ)=E[x1​(t)⋅x2​(t)]

由于数据是有限长，所以上面的信号的互相关可以使用下面公式进行估计：R^x1x2(τ)=1T−τ∫τTx1(t)⋅x2(t−τ)dτ\hat R_{x_1 x_2 } \left( \tau \right) = {1 \over {T - \tau }}\int_\tau ^T {x_1 \left( t \right) \cdot x_2 \left( {t - \tau } \right)d\tau }R^x1​x2​​(τ)=T−τ1​∫τT​x1​(t)⋅x2​(t−τ)dτ

其中TTT是观测时间长度。

■ 广义互相关

为了提高时延DDD的精度，通常需要对信号x1(t),x2(t)x_1 \left( t \right),x_2 \left( t \right)x1​(t),x2​(t)进行预先滤波。

下图中，对于信号x1(t),x2(t)x_1 \left( t \right),x_2 \left( t \right)x1​(t),x2​(t)通过各自的滤波器H1,H2H_1 ,H_2H1​,H2​进行滤波，获得y2(t),y2(t)y_2 \left( t \right),y_2 \left( t \right)y2​(t),y2​(t)。然后在计算它们的互相关，获得时延DDD的估计值。

▲ 对接收到的信号波形滤波、延迟、相乘、积分获得延迟参数

通过选择不同的滤波器H1(f),H2(f)H_1 \left( f \right),H_2 \left( f \right)H1​(f),H2​(f)，可以改善时延估计精度，通过滤波后数据的互相关峰值去逼近真实的时延DDD。

通过傅里叶变换可以获得接受信号的互功率谱密度函数Gx1x2(f)G_{x_1 x_2 } \left( f \right)Gx1​x2​​(f)：Rx1x2(τ)=∫−∞∞Gx2x2(f)ej2πfτdfR_{x_1 x_2 } \left( \tau \right) = \int_{ - \infty }^\infty {G_{x_2 x_2 } \left( f \right)e^{j2\pi f\tau } df}Rx1​x2​​(τ)=∫−∞∞​Gx2​x2​​(f)ej2πfτdf

x1(t),x2(t)x_1 \left( t \right),x_2 \left( t \right)x1​(t),x2​(t)之间的广义互相关可以通过下面公式给出：Ry1y2(g)(τ)=∫−∞∞ψ(f)Gx1x2(f)ej2πfτdfR_{y_1 y_2 }^{\left( g \right)} \left( \tau \right) = \int_{ - \infty }^\infty \psi \left( f \right)G_{x_1 x_2 } \left( f \right)e^{j2\pi f\tau } dfRy1​y2​(g)​(τ)=∫−∞∞​ψ(f)Gx1​x2​​(f)ej2πfτdf

其中：ψg(f)=H1(f)⋅H2∗(f)\psi _g \left( f \right) = H_1 \left( f \right) \cdot H_2^* \left( f \right)ψg​(f)=H1​(f)⋅H2∗​(f)

在实际应用中，信号的互功率谱Gx1x2(f)G_{x_1 x_2 } \left( f \right)Gx1​x2​​(f)只能通过观察的数据x1(t),x2(t)x_1 \left( t \right),x_2 \left( t \right)x1​(t),x2​(t)进行估算，所以：R^y1y2(g)(τ)=∫−∞∞ψ(f)G^x1x2(f)ej2πfτdf\hat R_{y_1 y_2 }^{\left( g \right)} \left( \tau \right) = \int_{ - \infty }^\infty {\psi \left( f \right)\hat G_{x_1 x_2 } \left( f \right)e^{j2\pi f\tau } df}R^y1​y2​(g)​(τ)=∫−∞∞​ψ(f)G^x1​x2​​(f)ej2πfτdf

■ 预处理

信号x1,x2x_1 ,x_2x1​,x2​之间的互相关：Rx1x2(τ)=a⋅Rs1s2(τ−D)+Rn1n2(τ)R_{x_1 x_2 } \left( \tau \right) = a \cdot R_{s_1 s_2 } \left( {\tau - D} \right) + R_{n_1 n_2 } \left( \tau \right)Rx1​x2​​(τ)=a⋅Rs1​s2​​(τ−D)+Rn1​n2​​(τ)

进行傅里叶变换，可得：Gx1x2(f)=a⋅Gs1s2(f)⋅e−j2πfD+Gn1n2(f)G_{x_1 x_2 } \left( f \right) = a \cdot G_{s_1 s_2 } \left( f \right) \cdot e^{ - j2\pi fD} + G_{n_1 n_2 } \left( f \right)Gx1​x2​​(f)=a⋅Gs1​s2​​(f)⋅e−j2πfD+Gn1​n2​​(f)

上面公众是的第一项，可以看做下面公式的傅里叶变换：Rx1x2(τ)=a⋅Rs1s21(τ)∗δ(t−D)R_{x1x2} \left( \tau \right) = a \cdot R_{s1s21} \left( \tau \right) * \delta \left( {t - D} \right)Rx1x2​(τ)=a⋅Rs1s21​(τ)∗δ(t−D)

对此可以解释为：信号的互相关是δ(t)\delta \left( t \right)δ(t)被信号的频谱进行扩展。

对于只有一个延迟问题并不严重。但对于多种信号的延迟，真实的互相关为：Rx1x2(τ)=Rs1s2(τ)∗∑iαiδ(τ−Di)R_{x1x2} \left( \tau \right) = R_{s1s2} \left( \tau \right) * \sum\limits_i^{} {\alpha _i \delta \left( {\tau - D_i } \right)}Rx1x2​(τ)=Rs1s2​(τ)∗i∑​αi​δ(τ−Di​)

1. Roth预处理

频谱的加权函数有Roth提出：ψR(f)=1Gx1x1(f)\psi _R \left( f \right) = {1 \over {G_{x1x1} \left( f \right)}}ψR​(f)=Gx1x1​(f)1​

这样互相关变成：

R^y1y2(R)(τ)=δ(τ−D)∗∫−∞∞α⋅Gs1s2(f)Gs1s2(f)+Gn1n2(f)⋅ej2πfτdf\hat R_{y1y2}^{\left( R \right)} \left( \tau \right) = \delta \left( {\tau - D} \right)*\int_{ - \infty }^\infty {{{\alpha \cdot G_{s1s2} \left( f \right)} \over {G_{s1s2} \left( f \right) + G_{n1n2} \left( f \right)}} \cdot e^{j2\pi f\tau } df}R^y1y2(R)​(τ)=δ(τ−D)∗∫−∞∞​Gs1s2​(f)+Gn1n2​(f)α⋅Gs1s2​(f)​⋅ej2πfτdf

当Gn1n2(f)G_{n1n2} \left( f \right)Gn1n2​(f)比较大的对应互相关进行了压制。

2.平滑向相关变换(SCOT)

加权函数为：KaTeX parse error: Double subscript at position 84: …) \cdot G_{x2} _̲{x2} \left( f \…

此时：R^y1y2(s)(τ)=∫−∞∞γ^x1x2(f)ej2πfτdf\hat R_{y1y2}^{\left( s \right)} \left( \tau \right) = \int_{ - \infty }^\infty {\hat \gamma _{x1x2} \left( f \right)e^{j2\pi f\tau } df}R^y1y2(s)​(τ)=∫−∞∞​γ^​x1x2​(f)ej2πfτdf

其中：KaTeX parse error: Undefined control sequence: \buildrel at position 38: …eft( f \right) \̲b̲u̲i̲l̲d̲r̲e̲l̲ ̲\Delta \over = …

当Gx1x1(f)=Gx2x2(f)G_{x1x1} \left( f \right) = G_{x2x2} \left( f \right)Gx1x1​(f)=Gx2x2​(f)，此时SCOT与Roth是等效的。

3.相位变换（PHAT）

PHAT使用加权函数：ψp(f)=1∣Gx1x2(f)∣\psi _p \left( f \right) = {1 \over {\left| {G_{x1x2} \left( f \right)} \right|}}ψp​(f)=∣Gx1x2​(f)∣1​

如果噪声n1(t),n2(t)n_1 \left( t \right),n_2 \left( t \right)n1​(t),n2​(t)不相关，即Gn1n2(f)=0G_{n1n2} \left( f \right) = 0Gn1n2​(f)=0，互相关为：Ry1y2(p)(τ)=δ(t−D)R_{y1y2}^{\left( p \right)} \left( \tau \right) = \delta \left( {t - D} \right)Ry1y2(p)​(τ)=δ(t−D)

其中：G^x1x2(f)∣Gx1x2(f)∣=ejθ(f)=ej2πfD{{\hat G_{x1x2} \left( f \right)} \over {\left| {G_{x1x2} \left( f \right)} \right|}} = e^{j\theta \left( f \right)} = e^{j2\pi fD}∣Gx1x2​(f)∣G^x1x2​(f)​=ejθ(f)=ej2πfD

当噪声之间不相关是，PHAT不会像其他预处理方法使得互相关扩散。

当G^x1x2(f)≠Gx1x2(f),θ(f)≠2πfD\hat G_{x1x2} \left( f \right) \ne G_{x1x2} \left( f \right),\theta \left( f \right) \ne 2\pi fDG^x1x2​(f)​=Gx1x2​(f),θ(f)​=2πfD

Ry1y2(p)(τ)R_{y1y2}^{\left( p \right)} \left( \tau \right)Ry1y2(p)​(τ)不等于delta函数。

4. Eckart 滤波

加权函数为：ψE(f)=αGs1s2(f)Gn1n1(f)⋅Gn2n2(f)\psi _E \left( f \right) = {{\alpha G_{s1s2} \left( f \right)} \over {G_{n1n1} \left( f \right) \cdot G_{n2n2} \left( f \right)}}ψE​(f)=Gn1n1​(f)⋅Gn2n2​(f)αGs1s2​(f)​

5. Hannon与Thomson

Ry1y2(HT)(τ)=∫−∞∞G^x1x2(f)⋅1∣Gx1x2(f)∣⋅∣γ12(f)∣2[1−∣γ12(f)∣2]ej2πfτdfR_{y1y2}^{\left( {HT} \right)} \left( \tau \right) = \int_{ - \infty }^\infty {\hat G_{x1x2} \left( f \right) \cdot {1 \over {\left| {G_{x1x2} \left( f \right)} \right|}} \cdot {{\left| {\gamma _{12} \left( f \right)} \right|^2 } \over {\left[ {1 - \left| {\gamma _{12} \left( f \right)} \right|^2 } \right]}}e^{j2\pi f\tau } df}Ry1y2(HT)​(τ)=∫−∞∞​G^x1x2​(f)⋅∣Gx1x2​(f)∣1​⋅[1−∣γ12​(f)∣2]∣γ12​(f)∣2​ej2πfτdf

对应的频谱加权函数为：

ψHT(f)=1∣Gx1x2(f)∣⋅∣γ12(f)∣2[1−∣γ12(f)∣2]\psi _{HT} \left( f \right) = {1 \over {\left| {G_{x1x2} \left( f \right)} \right|}} \cdot {{\left| {\gamma _{12} \left( f \right)} \right|^2 } \over {\left[ {1 - \left| {\gamma _{12} \left( f \right)} \right|^2 } \right]}}ψHT​(f)=∣Gx1x2​(f)∣1​⋅[1−∣γ12​(f)∣2]∣γ12​(f)∣2​

■ 对于低信噪比(SNR)下的最大似然估计

在信噪比低的情况下：Gs1s1(f)Gn1n1(f)≪1,Gs1s1(f)Gn2n2≪1{{G_{s1s1} \left( f \right)} \over {G_{n1n1} \left( f \right)}} \ll 1,\,\,\,\,{{G_{s1s1} \left( f \right)} \over {G_{n2n2} }} \ll 1Gn1n1​(f)Gs1s1​(f)​≪1,Gn2n2​Gs1s1​(f)​≪1

ψHT(f)≅Gs1s1(f)Gn1n1(f)⋅Gn2n2(f)=ψE(f)\psi _{HT} \left( f \right) \cong {{G_{s1s1} \left( f \right)} \over {G_{n1n1} \left( f \right) \cdot G_{n2n2} \left( f \right)}} = \psi _E \left( f \right)ψHT​(f)≅Gn1n1​(f)⋅Gn2n2​(f)Gs1s1​(f)​=ψE​(f)

如果：Gn1n1(f)=Gn2n2(f)=Gnn(f)G_{n1n1} \left( f \right) = G_{n2n2} \left( f \right) = G_{nn} \left( f \right)Gn1n1​(f)=Gn2n2​(f)=Gnn​(f)

那么：ψHT(f)≅Gs1s1(f)Gnn(f)[ψs(f)]=[Gs1s1(f)Gnn(f)]2ψp(f)\psi _{HT} \left( f \right) \cong {{G_{s1s1} \left( f \right)} \over {G_{nn} \left( f \right)}}\left[ {\psi _s \left( f \right)} \right] = \left[ {{{G_{s1s1} \left( f \right)} \over {G_{nn} \left( f \right)}}} \right]^2 \psi _p \left( f \right)ψHT​(f)≅Gnn​(f)Gs1s1​(f)​[ψs​(f)]=[Gnn​(f)Gs1s1​(f)​]2ψp​(f)

■ 结论

HT预处理适应于通常情况下的延迟估计；在SNR低的情况下，HT与Eckart，互相关等效；如果互相关随着时间缓慢变化，最大似然（ML）估计中的滤波器也是随着时间变化的。