python - 如何缩放基于 FFT 的互相关，使其峰值等于 Pearson 的 rho

Question

问题描述

FFT 可用于计算两个信号或图像之间的互相关。要确定两个信号 A 和 B 之间的延迟或滞后，只需找到以下峰值:
IFFT(FFT(A)*共轭(FFT(B)))


但是，峰值的幅度与各个信号的频谱幅度有关。从而确定Pearson correlation (rho) ，这个峰值的幅度必须由两个信号中的总能量来衡量。

执行此操作的一种方法是 normalize by the geometric mean of the individual autocorrelations .这给出了 rho 的合理近似值，尤其是当样本之间的延迟很小但不是精确值时。

我认为出现此错误的原因是 PIL 逊相关性仅针对信号的重叠部分定义，而归一化因子(两个自相关峰的几何平均值)包括非重叠部分的贡献。我考虑了两种解决此问题并通过 FFT 为 rho 生成精确值的方法。在第一个(下面称为 rho_exact_1)中，我将样本修剪到它们的重叠部分，并从中计算归一化因子。在第二个(下面称为 rho_exact_2)中，我计算了信号重叠部分中包含的测量值的分数，并将完全自相关归一化因子乘以该分数。

都不行!下图显示了使用基于 DFT 的互相关计算 Pearson 的 rho 的三种方法的图。仅显示了互相关峰的区域。每个估计都接近正确值 1.0，但不等于它。

我用来执行计算的代码如下。我使用一个简单的正弦波作为示例信号。我注意到，如果我使用方波(占空比不一定为 50%)，方法的误差会发生变化。

谁能解释一下这是怎么回事？

一个工作示例

import numpy as np
from matplotlib import pyplot as plt

# make a time vector w/ 256 points
# and a source signal
N_cycles = 10.0
N_points = 256.0

t = np.arange(0,N_cycles*np.pi,np.pi*N_cycles/N_points)
signal = np.sin(t)

use_rect = False
if use_rect:
    threshold = -0.75
    signal[np.where(signal>=threshold)]=1.0
    signal[np.where(signal<threshold)]=-1.0

# normalize the signal (not technically
# necessary for this example, but required
# for measuring correlation of physically
# different signals)
signal = signal/signal.std()

# generate two samples of the signal
# with a temporal offset:
N = 128
offset = 5
sample_1 = signal[:N]
sample_2 = signal[offset:N+offset]

# determine the offset through cross-
# correlation
xc_num = np.abs(np.fft.ifft(np.fft.fft(sample_1)*np.fft.fft(sample_2).conjugate()))
offset_estimate = np.argmax(xc_num)
if offset_estimate>N//2:
    offset_estimate = offset_estimate - N


# for an approximate estimate of Pearson's
# correlation, we normalize by the RMS
# of individual autocorrelations:
autocorrelation_1 = np.abs(np.fft.ifft(np.fft.fft(sample_1)*np.fft.fft(sample_1).conjugate()))
autocorrelation_2 = np.abs(np.fft.ifft(np.fft.fft(sample_2)*np.fft.fft(sample_2).conjugate()))
xc_denom_approx = np.sqrt(np.max(autocorrelation_1))*np.sqrt(np.max(autocorrelation_2))
rho_approx = xc_num/xc_denom_approx
print 'rho_approx',np.max(rho_approx)

# this is an approximation because we've
# included autocorrelation of the whole samples
# instead of just the overlapping portion;
# using cropped versions of the samples should
# yield the correct correlation:
sample_1_cropped = sample_1[offset:]
sample_2_cropped = sample_2[:-offset]

# these should be identical vectors:
assert np.all(sample_1_cropped==sample_2_cropped)

# compute autocorrelations of cropped samples
# and corresponding value for rho
autocorrelation_1_cropped = np.abs(np.fft.ifft(np.fft.fft(sample_1_cropped)*np.fft.fft(sample_1_cropped).conjugate()))
autocorrelation_2_cropped = np.abs(np.fft.ifft(np.fft.fft(sample_2_cropped)*np.fft.fft(sample_2_cropped).conjugate()))
xc_denom_exact_1 = np.sqrt(np.max(autocorrelation_1_cropped))*np.sqrt(np.max(autocorrelation_2_cropped))
rho_exact_1 = xc_num/xc_denom_exact_1
print 'rho_exact_1',np.max(rho_exact_1)

# alternatively we could try to use the
# whole sample autocorrelations and just
# scale by the number of pixels used to
# compute the numerator:
scaling_factor = float(len(sample_1_cropped))/float(len(sample_1))
rho_exact_2 = xc_num/(xc_denom_approx*scaling_factor)
print 'rho_exact_2',np.max(rho_exact_2)

# finally a sanity check: is rho actually 1.0
# for the two signals:
rho_corrcoef = np.corrcoef(sample_1_cropped,sample_2_cropped)[0,1]
print 'rho_corrcoef',rho_corrcoef

x = np.arange(len(rho_approx))
plt.plot(x,rho_approx,label='FFT rho_approx')
plt.plot(x,rho_exact_1,label='FFT rho_exact_1')
plt.plot(x,rho_exact_2,label='FFT rho_exact_2')
plt.plot(x,np.ones(len(x))*rho_corrcoef,'k--',label='Pearson rho')
plt.legend()
plt.ylim((.75,1.25))
plt.xlim((0,20))
plt.show()

Answer 1

两个N周期离散信号F和G之间的归一化互相关定义为:

由于分子是两个向量(F 和 G_x)之间的点积，而分母是这两个向量范数的乘积，标量 r_x 必须确实位于 -1 和 +1 之间，它是余弦矢量之间的角度(参见 there )。如果向量F和G_x是对齐的，那么r_x=1。如果 r_x=1，则向量 F 和 G_x 由于三角不等式而对齐。为确保这些属性，分子处的向量必须与分母处的向量匹配。

可以使用离散傅立叶变换一次计算所有分子。实际上，该变换将卷积转换为傅立叶空间中的逐点乘积。这就是为什么在您执行的测试中不同的估计归一化互相关不是 1。

1. 对于第一个测试“approx”，sample_1 和sample_2 都是从周期信号中提取的。两者的长度相同，但长度不是周期的倍数，因为它是 2.5 个周期 (5pi)(下图)。结果，由于 dft 执行相关，就好像它们在周期性信号中一样，发现 sample_1 和 sample_2 不是完全相关的并且 r_x<1。

2. 对于第二个测试 rho_exact_1，卷积是在长度为 N=128 的信号上执行的，但是分母处的范数是在大小为 N-offset=128- 的截断向量上计算的5.结果，r_x 的属性丢失了。此外，必须注意所提出的卷积和范数未归一化:计算出的范数和卷积乘积与所考虑向量的点数全局成正比。因此，与前一种情况相比，截断向量的范数略低，并且 r_x 增加:随着 offset 增加，可能会遇到大于 1 的值。

3. 对于第三个测试rho_exact_2，引入了一个缩放因子来尝试纠正第一个测试:r_x 的属性也丢失了，并且在缩放时可能会遇到大于 1 的值因子大于一。

尽管如此，numpy 的函数 corrcoef() 实际上为截断信号计算了一个等于 1 的 r_x。事实上，这些信号完全相同!使用 DFT 可以获得相同的结果:

xc_num_cropped = np.abs(np.fft.ifft(np.fft.fft(sample_1_cropped)*np.fft.fft(sample_2_cropped).conjugate()))
autocorrelation_1_cropped = np.abs(np.fft.ifft(np.fft.fft(sample_1_cropped)*np.fft.fft(sample_1_cropped).conjugate()))
autocorrelation_2_cropped = np.abs(np.fft.ifft(np.fft.fft(sample_2_cropped)*np.fft.fft(sample_2_cropped).conjugate()))
xc_denom_exact_11 = np.sqrt(np.max(autocorrelation_1_cropped))*np.sqrt(np.max(autocorrelation_2_cropped))
rho_exact_11 = xc_num_cropped/xc_denom_exact_11
print 'rho_exact_11',np.max(rho_exact_11)  

要为用户提供重要的 r_x 值，您可以坚持第一个测试提供的值，如果帧的长度不是周期的倍数，则对于相同的周期信号，该值可能低于 1。为了纠正这个缺点，估计的偏移量也可以被检索并用于构建两个相同长度的裁剪信号。必须重新运行整个相关过程以获得 r_x 的新值，这不会因为裁剪帧的长度不是周期的倍数而受到影响。

最后，如果 DFT 是一次计算所有 x 值的分子卷积的非常有效的方法，则可以使用 numpy.linalg.norm 将分母有效地计算为向量的 2 范数。由于如果第一次相关成功，裁剪信号的 argmax(r_x) 可能为零，因此使用点积 `sample_1_cropped.dot(sample_2_cropped) 计算 r_0 就足够了。