Автокорреляция многомерного массива в NumPy

Question

У меня есть двумерный массив, то есть массив последовательностей, которые также являются массивами. Для каждой последовательности я хотел бы рассчитать автокорреляцию, чтобы для массива (5,4) я получил 5 результатов или массив измерения (5,7).

Я знаю, что мог бы просто пройтись по первому измерению, но это медленно и в крайнем случае. Есть ли другой способ?

Спасибо!

EDIT:

На основании выбранного ответа плюс комментарий от mtrw у меня есть следующая функция:

def xcorr(x):
  """FFT based autocorrelation function, which is faster than numpy.correlate"""
  # x is supposed to be an array of sequences, of shape (totalelements, length)
  fftx = fft(x, n=(length*2-1), axis=1)
  ret = ifft(fftx * np.conjugate(fftx), axis=1)
  ret = fftshift(ret, axes=1)
  return ret

Обратите внимание, что длина - это глобальная переменная в моем коде, поэтому обязательно объявите ее. Я также не ограничивал результат действительными числами, поскольку мне также необходимо учитывать комплексные числа.

Answer 1

Использование автокорреляции на основе БПФ :

import numpy
from numpy.fft import fft, ifft

data = numpy.arange(5*4).reshape(5, 4)
print data
##[[ 0  1  2  3]
## [ 4  5  6  7]
## [ 8  9 10 11]
## [12 13 14 15]
## [16 17 18 19]]
dataFT = fft(data, axis=1)
dataAC = ifft(dataFT * numpy.conjugate(dataFT), axis=1).real
print dataAC
##[[   14.     8.     6.     8.]
## [  126.   120.   118.   120.]
## [  366.   360.   358.   360.]
## [  734.   728.   726.   728.]
## [ 1230.  1224.  1222.  1224.]]

Меня немного смущает ваше утверждение о том, что ответ имеет размерность (5, 7), поэтому, возможно, есть что-то важное, чего я не понимаю.

РЕДАКТИРОВАТЬ: По предложению mtrw, мягкая версия, которая не распространяется:

import numpy
from numpy.fft import fft, ifft

data = numpy.arange(5*4).reshape(5, 4)
padding = numpy.zeros((5, 3))
dataPadded = numpy.concatenate((data, padding), axis=1)
print dataPadded
##[[  0.   1.   2.   3.   0.   0.   0.   0.]
## [  4.   5.   6.   7.   0.   0.   0.   0.]
## [  8.   9.  10.  11.   0.   0.   0.   0.]
## [ 12.  13.  14.  15.   0.   0.   0.   0.]
## [ 16.  17.  18.  19.   0.   0.   0.   0.]]
dataFT = fft(dataPadded, axis=1)
dataAC = ifft(dataFT * numpy.conjugate(dataFT), axis=1).real
print numpy.round(dataAC, 10)[:, :4]
##[[   14.     8.     3.     0.     0.     3.     8.]
## [  126.    92.    59.    28.    28.    59.    92.]
## [  366.   272.   179.    88.    88.   179.   272.]
## [  734.   548.   363.   180.   180.   363.   548.]
## [ 1230.   920.   611.   304.   304.   611.   920.]]

Должен быть более эффективный способ сделать это, особенно потому, что автокорреляция симметрична, и я этим не пользуюсь.

Answer 2

Для действительно больших массивов становится важным иметь n = 2 ** p, где p - целое число.Это сэкономит вам огромное количество времени.Например:

def xcorr(x):
  l = 2 ** int(np.log2(length * 2 - 1))
  fftx = fft(x, n = l, axis = 1)
  ret = ifft(fftx * np.conjugate(fftx), axis = 1)
  ret = fftshift(ret, axes=1)
  return ret

Это может привести к ошибкам при переходе.Однако для больших массивов автокорреляция должна быть незначительной вблизи краев.

Answer 3

Может быть, это просто предпочтение, но я хотел следовать из определения. Лично мне так легче следовать. Это моя реализация для произвольного массива nd.

from itertools import product
from numpy import empty, roll</p>

<code>def autocorrelate(x):
    """
    Compute the multidimensional autocorrelation of an nd array.
    input: an nd array of floats
    output: an nd array of autocorrelations
    """

    # used for transposes
    t = roll(range(x.ndim), 1)

    # pairs of indexes
    # the first is for the autocorrelation array
    # the second is the shift
    ii = [list(enumerate(range(1, s - 1))) for s in x.shape]

    # initialize the resulting autocorrelation array
    acor = empty(shape=[len(s0) for s0 in ii])

    # iterate over all combinations of directional shifts
    for i in product(*ii):
        # extract the indexes for
        # the autocorrelation array 
        # and original array respectively
        i1, i2 = asarray(i).T

        x1 = x.copy()
        x2 = x.copy()

        for i0 in i2:
            # clip the unshifted array at the end
            x1 = x1[:-i0]
            # and the shifted array at the beginning
            x2 = x2[i0:]

            # prepare to do the same for 
            # the next axis
            x1 = x1.transpose(t)
            x2 = x2.transpose(t)

        # normalize shifted and unshifted arrays
        x1 -= x1.mean()
        x1 /= x1.std()
        x2 -= x2.mean()
        x2 /= x2.std()

        # compute the autocorrelation directly
        # from the definition
        acor[tuple(i1)] = (x1 * x2).mean()

    return acor
</code>