问题描述
很可能是我现在让自己看起来很愚蠢，但以下是从基于 fft 的自相关获得间隔 [1,-1] 缩放输出的正确方法吗?缩放应该是 numpy.correlate(x,x, mode="same") 将结果缩放到 [1, -1] 区间.
It could well be i am making myself look very stupid just now but is the following the correct way to get interval[1,-1] scaled output from fft based autocorrelation? The scaling should be what numpy.correlate(x,x, mode="same") does to scale the results to a [1, -1] interval.
def autocorrelate(x):
fftx = fft(x)
fftx_mean = np.mean(fftx)
fftx_std = np.std(fftx)
ffty = np.conjugate(fftx)
ffty_mean = np.mean(ffty)
ffty_std = np.std(ffty)
result = ifft((fftx - fftx_mean) * (ffty - ffty_mean))
result = fftshift(result)
return [i/(fftx_std * ffty_std) for i in result.real]
我已经运行了一些测试数据，它确实看起来像它应该做的那样，但我不确定我没有搞砸一些事情，只是不小心得到了一些正确的结果;)
i have run some test data and it sure looks like it does what it should but i am not perfectly sure i haven't screwed something up and just accidentally get somewhat correct results ;)
推荐答案
Maple 的 AutoCorrelation 函数好像使用了定义
Maple's AutoCorrelation function seems to be using the definition
def AutoCorrelation(x):
x = np.asarray(x)
y = x-x.mean()
result = np.correlate(y, y, mode='full')
result = result[len(result)//2:]
result /= result[0]
return result
In [189]: AutoCorrelation([1,2,1,2,1,2,1,2])
Out[189]: array([ 1. , -0.875, 0.75 , -0.625, 0.5 , -0.375, 0.25 , -0.125])
现在，看看我们是否可以使用 FFT 重现这个结果会很有趣.NumPy 的 np.fft.fft 是一个 周期性卷积，而 np.correlate 是一个线性卷积.要使用 np.fft.fft，我们需要添加足够的零填充以使计算基本线性:
Now, it would be interesting to see if we can reproduce this result using FFT. NumPy's np.fft.fft is a periodic convolution, while np.correlate is a linear convolution. To use np.fft.fft, we need to add enough zero-padding to make the calculation essentially linear:
def autocorrelation(x):
"""
Compute autocorrelation using FFT
"""
x = np.asarray(x)
N = len(x)
x = x-x.mean()
s = fft.fft(x, N*2-1)
result = np.real(fft.ifft(s * np.conjugate(s), N*2-1))
result = result[:N]
result /= result[0]
return result
这里有一些测试确认 AutoCorrelation 和 autocorrelation 同意并返回与 Maple 的 AutoCorrelation 函数返回的值相同的值——至少对于我知道的有限示例约.
Here are some tests which confirm that AutoCorrelation and autocorrelation agree and return the same values as those returned by Maple's AutoCorrelation function -- at least for the limited examples I know about.
import numpy as np
fft = np.fft
def autocorrelation(x):
"""
Compute autocorrelation using FFT
The idea comes from
http://dsp.stackexchange.com/a/1923/4363 (Hilmar)
"""
x = np.asarray(x)
N = len(x)
x = x-x.mean()
s = fft.fft(x, N*2-1)
result = np.real(fft.ifft(s * np.conjugate(s), N*2-1))
result = result[:N]
result /= result[0]
return result
def AutoCorrelation(x):
x = np.asarray(x)
y = x-x.mean()
result = np.correlate(y, y, mode='full')
result = result[len(result)//2:]
result /= result[0]
return result
def autocorrelate(x):
fftx = fft.fft(x)
fftx_mean = np.mean(fftx)
fftx_std = np.std(fftx)
ffty = np.conjugate(fftx)
ffty_mean = np.mean(ffty)
ffty_std = np.std(ffty)
result = fft.ifft((fftx - fftx_mean) * (ffty - ffty_mean))
result = fft.fftshift(result)
return [i / (fftx_std * ffty_std) for i in result.real]
np.set_printoptions(precision=3, suppress=True)
"""
These tests come from
http://www.maplesoft.com/support/help/Maple/view.aspx?path=Statistics/AutoCorrelation
http://www.maplesoft.com/support/help/Maple/view.aspx?path=updates%2fMaple15%2fcomputation
"""
tests = [
([1,2,1,2,1,2,1,2], [1,-0.875,0.75,-0.625,0.5,-0.375,0.25,-0.125]),
([1,-1,1,-1], [1, -0.75, 0.5, -0.25]),
]
for x, answer in tests:
x = np.array(x)
answer = np.array(answer)
# print(autocorrelate(x))
print(autocorrelation(x))
print(AutoCorrelation(x))
assert np.allclose(AutoCorrelation(x), answer)
print
"""
Test that autocorrelation() agrees with AutoCorrelation()
"""
for i in range(1000):
x = np.random.random(np.random.randint(2,100))*100
assert np.allclose(autocorrelation(x), AutoCorrelation(x))
这篇关于使用 numpy.fft 对自相关进行统计缩放的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！