3

很可能我刚才让自己看起来很愚蠢，但以下是从基于 fft 的自相关获得间隔 [1，-1] 缩放输出的正确方法吗？缩放应该是 numpy.correlate(x,x, mode="same") 所做的，以将结果缩放到 [1, -1] 间隔。

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]

我已经运行了一些测试数据，它确实看起来像它应该做的，但我不完全确定我没有搞砸什么，只是不小心得到了一些正确的结果；）

4

1 回答 1

7

Maple 的 AutoCorrelation 函数似乎使用了定义

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 重现这个结果会很有趣。NumPynp.fft.fft是周期性卷积，np.correlate而是线性卷积。要使用np.fft.fft，我们需要添加足够的零填充以使计算基本上是线性的：

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 函数返回的值相同的值——至少对于我所知道的有限示例而言。

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))
于 2013-04-16T19:10:12.740 回答