使用numpy.fft进行自相关的统计缩放

Question

使用numpy.fft进行自相关的统计缩放

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]

我已经运行了一些测试数据，它确实看起来像是在做它应该做的事情，但我不完全确定我没有搞砸，并且只是偶然得到了一些正确的结果 ;)

- Skazarok

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- unutbu · Accepted Answer

据我了解，Maple的自相关函数似乎使用的是定义式。

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，我们需要添加足够的零填充来使计算基本上呈线性：

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))