问题描述

FFT可用于计算两个信号或图像之间的互相关。确定两个信号之间的延迟或滞后 A. 和 ,它足以定位以下峰值: IFFT(FFT(A)*conjugate(FFT(B)))

Pearson correlation (rho) ,该峰值的幅度必须根据两个信号中的总能量进行缩放。

一种方法是 normalize by the geometric mean of the individual autocorrelations rho公司

我认为这个错误的原因是皮尔逊相关性仅定义为信号的重叠部分,而归一化因子(两个自相关峰的几何平均值)包括非重叠部分的贡献。我考虑了两种方法来解决这个问题,并为 rho公司 rho_exact_1 下面),我将样本裁剪到其重叠部分,并根据这些部分计算归一化因子。在第二个(调用 rho_exact_2 下面),我计算了信号重叠部分中包含的测量分数,并将全自相关归一化因子乘以该分数。

两者都不起作用!下图显示了计算皮尔逊指数的三种方法的曲线图 rho公司

我用于执行计算的代码如下。我使用了一个简单的正弦波作为示例信号。我注意到,如果我使用方波(占空比不一定为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()