使用 numpy.fft 对自相关进行统计缩放 [英] Statistical Scaling of autocorrelation using numpy.fft

问题描述

很可能是我现在让自己看起来很愚蠢，但以下是从基于 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 对自相关进行统计缩放的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！