广义互相关的公式,这一文都搜集全了

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

 

■ 互相关估计时延


在空间上相距一定距离 L L L的两个声音传感器获取远方声源发出的声音信号可以表示为: x 1 ( t ) = s 1 ( t ) + n 1 ( t ) x_1 \left( t \right) = s_1 \left( t \right) + n_1 \left( t \right) x1(t)=s1(t)+n1(t) x 2 ( t ) = a ⋅ s 1 ( t + D ) + n 2 ( t ) x_2 \left( t \right) = a \cdot s_1 \left( {t + D} \right) + n_2 \left( t \right) x2(t)=as1(t+D)+n2(t)
其中 s 1 ( t ) , n 1 ( t ) , n 2 ( t ) s_1 \left( t \right),n_1 \left( t \right),n_2 \left( t \right) s1(t),n1(t),n2(t)是实信号,可以看成联合平稳随机过程。信号 s 1 ( t ) s_1 \left( t \right) s1(t)与两个噪声 n 1 ( t ) , n 2 ( t ) n_1 \left( t \right),n_2 \left( t \right) n1(t),n2(t)之间是不相关的。

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

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

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

一种常用的计算时间延迟 D D D的方法是通过计算两个传感器信号的互相关来获得的。基于时间延迟 D D D,相对于传感器的信号波达方向也可以计算出来: R x 1 , x 2 ( τ ) = E [ x 1 ( t ) ⋅ x 2 ( 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 ^ x 1 x 2 ( τ ) = 1 T − τ ∫ τ T x 1 ( t ) ⋅ x 2 ( 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^x1x2(τ)=Tτ1τTx1(t)x2(tτ)dτ
其中 T T T是观测时间长度。

 

■ 广义互相关


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

下图中,对于信号 x 1 ( t ) , x 2 ( t ) x_1 \left( t \right),x_2 \left( t \right) x1(t),x2(t)通过各自的滤波器 H 1 , H 2 H_1 ,H_2 H1,H2进行滤波,获得 y 2 ( t ) , y 2 ( t ) y_2 \left( t \right),y_2 \left( t \right) y2(t),y2(t)。然后在计算它们的互相关,获得时延 D D D的估计值。

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

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

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

通过傅里叶变换可以获得接受信号的互功率谱密度函数 G x 1 x 2 ( f ) G_{x_1 x_2 } \left( f \right) Gx1x2(f) R x 1 x 2 ( τ ) = ∫ − ∞ ∞ G x 2 x 2 ( f ) e j 2 π f τ d f R_{x_1 x_2 } \left( \tau \right) = \int_{ - \infty }^\infty {G_{x_2 x_2 } \left( f \right)e^{j2\pi f\tau } df} Rx1x2(τ)=Gx2x2(f)ej2πfτdf

x 1 ( t ) , x 2 ( t ) x_1 \left( t \right),x_2 \left( t \right) x1(t),x2(t)之间的广义互相关可以通过下面公式给出: R y 1 y 2 ( g ) ( τ ) = ∫ − ∞ ∞ ψ ( f ) G x 1 x 2 ( f ) e j 2 π f τ d f R_{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 } df Ry1y2(g)(τ)=ψ(f)Gx1x2(f)ej2πfτdf
其中: ψ g ( f ) = H 1 ( f ) ⋅ H 2 ∗ ( f ) \psi _g \left( f \right) = H_1 \left( f \right) \cdot H_2^* \left( f \right) ψg(f)=H1(f)H2(f)
在实际应用中,信号的互功率谱 G x 1 x 2 ( f ) G_{x_1 x_2 } \left( f \right) Gx1x2(f)只能通过观察的数据 x 1 ( t ) , x 2 ( t ) x_1 \left( t \right),x_2 \left( t \right) x1(t),x2(t)进行估算,所以: R ^ y 1 y 2 ( g ) ( τ ) = ∫ − ∞ ∞ ψ ( f ) G ^ x 1 x 2 ( f ) e j 2 π f τ d f \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^y1y2(g)(τ)=ψ(f)G^x1x2(f)ej2πfτdf

 

■ 预处理


信号 x 1 , x 2 x_1 ,x_2 x1,x2之间的互相关: R x 1 x 2 ( τ ) = a ⋅ R s 1 s 2 ( τ − D ) + R n 1 n 2 ( τ ) 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) Rx1x2(τ)=aRs1s2(τD)+Rn1n2(τ)
进行傅里叶变换,可得: G x 1 x 2 ( f ) = a ⋅ G s 1 s 2 ( f ) ⋅ e − j 2 π f D + G n 1 n 2 ( 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) Gx1x2(f)=aGs1s2(f)ej2πfD+Gn1n2(f)
上面公众是的第一项,可以看做下面公式的傅里叶变换: R x 1 x 2 ( τ ) = a ⋅ R s 1 s 21 ( τ ) ∗ δ ( t − D ) R_{x1x2} \left( \tau \right) = a \cdot R_{s1s21} \left( \tau \right) * \delta \left( {t - D} \right) Rx1x2(τ)=aRs1s21(τ)δ(tD)
对此可以解释为:信号的互相关是 δ ( t ) \delta \left( t \right) δ(t)被信号的频谱进行扩展。

对于只有一个延迟问题并不严重。但对于多种信号的延迟,真实的互相关为: R x 1 x 2 ( τ ) = R s 1 s 2 ( τ ) ∗ ∑ i α i δ ( τ − D i ) 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 ) = 1 G x 1 x 1 ( f ) \psi _R \left( f \right) = {1 \over {G_{x1x1} \left( f \right)}} ψR(f)=Gx1x1(f)1
这样互相关变成:
R ^ y 1 y 2 ( R ) ( τ ) = δ ( τ − D ) ∗ ∫ − ∞ ∞ α ⋅ G s 1 s 2 ( f ) G s 1 s 2 ( f ) + G n 1 n 2 ( f ) ⋅ e j 2 π f τ d f \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
G n 1 n 2 ( 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 ^ y 1 y 2 ( s ) ( τ ) = ∫ − ∞ ∞ γ ^ x 1 x 2 ( f ) e j 2 π f τ d f \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 = …
G x 1 x 1 ( f ) = G x 2 x 2 ( f ) G_{x1x1} \left( f \right) = G_{x2x2} \left( f \right) Gx1x1(f)=Gx2x2(f),此时SCOT与Roth是等效的。

 

3.相位变换(PHAT)

PHAT使用加权函数: ψ p ( f ) = 1 ∣ G x 1 x 2 ( f ) ∣ \psi _p \left( f \right) = {1 \over {\left| {G_{x1x2} \left( f \right)} \right|}} ψp(f)=Gx1x2(f)1
如果噪声 n 1 ( t ) , n 2 ( t ) n_1 \left( t \right),n_2 \left( t \right) n1(t),n2(t)不相关,即 G n 1 n 2 ( f ) = 0 G_{n1n2} \left( f \right) = 0 Gn1n2(f)=0,互相关为: R y 1 y 2 ( p ) ( τ ) = δ ( t − D ) R_{y1y2}^{\left( p \right)} \left( \tau \right) = \delta \left( {t - D} \right) Ry1y2(p)(τ)=δ(tD)
其中: G ^ x 1 x 2 ( f ) ∣ G x 1 x 2 ( f ) ∣ = e j θ ( f ) = e j 2 π f D {{\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 ^ x 1 x 2 ( f ) ≠ G x 1 x 2 ( f ) , θ ( f ) ≠ 2 π f D \hat G_{x1x2} \left( f \right) \ne G_{x1x2} \left( f \right),\theta \left( f \right) \ne 2\pi fD G^x1x2(f)=Gx1x2(f),θ(f)=2πfD
R y 1 y 2 ( p ) ( τ ) R_{y1y2}^{\left( p \right)} \left( \tau \right) Ry1y2(p)(τ)不等于delta函数。
 

4. Eckart 滤波

加权函数为: ψ E ( f ) = α G s 1 s 2 ( f ) G n 1 n 1 ( f ) ⋅ G n 2 n 2 ( 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

R y 1 y 2 ( H T ) ( τ ) = ∫ − ∞ ∞ G ^ x 1 x 2 ( f ) ⋅ 1 ∣ G x 1 x 2 ( f ) ∣ ⋅ ∣ γ 12 ( f ) ∣ 2 [ 1 − ∣ γ 12 ( f ) ∣ 2 ] e j 2 π f τ d f R_{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)2ej2πfτdf
对应的频谱加权函数为:
ψ H T ( f ) = 1 ∣ G x 1 x 2 ( 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)下的最大似然估计


在信噪比低的情况下: G s 1 s 1 ( f ) G n 1 n 1 ( f ) ≪ 1 ,      G s 1 s 1 ( f ) G n 2 n 2 ≪ 1 {{G_{s1s1} \left( f \right)} \over {G_{n1n1} \left( f \right)}} \ll 1,\,\,\,\,{{G_{s1s1} \left( f \right)} \over {G_{n2n2} }} \ll 1 Gn1n1(f)Gs1s1(f)1,Gn2n2Gs1s1(f)1
ψ H T ( f ) ≅ G s 1 s 1 ( f ) G n 1 n 1 ( f ) ⋅ G n 2 n 2 ( 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)

如果: G n 1 n 1 ( f ) = G n 2 n 2 ( f ) = G n n ( f ) G_{n1n1} \left( f \right) = G_{n2n2} \left( f \right) = G_{nn} \left( f \right) Gn1n1(f)=Gn2n2(f)=Gnn(f)
那么: ψ H T ( f ) ≅ G s 1 s 1 ( f ) G n n ( f ) [ ψ s ( f ) ] = [ G s 1 s 1 ( f ) G n n ( 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)估计中的滤波器也是随着时间变化的。
©️2020 CSDN 皮肤主题: Age of Ai 设计师:meimeiellie 返回首页