Статистическое масштабирование автокорреляции с использованием numpy.fft

Question

Может быть, я делаю себя очень глупым сейчас, но это правильный способ получить интервал [1, -1] масштабированный вывод из автокорреляции на основе fft? Масштабирование должно быть таким, как 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] 

я запустить некоторые тестовые данные, и он уверен, выглядит, как он делает то, что он должен, но я не совсем уверен, что я не ввинчивается что-то и просто случайно получить несколько правильных результатов;)

Answer 1

Мэйпл-х автокорреляционная функция, как представляется, используя определение

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

Теперь, было бы интересно посмотреть, если мы можем воспроизвести этот результат с помощью быстрого преобразования Фурье. NumPy's np.fft.fft - это periodic convolution, а 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))