Блог

Главная - Вопросы по программированию - Автокорреляция нескольких временных рядов в Matlab с использованием FFT

Автокорреляция нескольких временных рядов в Matlab с использованием FFT

#matlab #time-series #fft #correlation

#matlab #временные ряды #fft #корреляция

Вопрос:

Функция Matlab autocorr (документация здесь ) вычисляет выборочную автокорреляцию одного временного ряда с использованием алгоритма быстрого преобразования Фурье (FFT):

 nFFT = 2^(nextpow2(length(y)) 1);
F = fft(y-mean(y),nFFT);
F = F.*conj(F);
acf = ifft(F);
acf = acf(1:(numLags 1)); % Retain non-negative lags
acf = acf./acf(1); % Normalize
acf = real(acf);

Предположим, у меня есть несколько реализаций одного и того же случайного процесса (временных рядов), и я хотел бы рассчитать выборочную автокорреляцию. Наивным подходом было бы вызывать autocorr для каждого временного ряда и усреднять корреляцию при каждом запаздывании. Однако лучше рассматривать все пары точек, разделенных запаздыванием k одновременно. Я реализовал это, построив матрицу (n-iLag)*nSamples на 2 из всех пар, разделенных размером задержки iLag , а затем вычислив корреляцию выборки. 

 % input: samplesMat (nSamples (times series) by n (time points) matrix)

sizeMat = size(samplesMat);
nSamples = sizeMat(1); 
n = sizeMat(2);

ACF = zeros(n,1);

muHat = mean(mean(samplesMat)); % sample mean

for iLag = 0:(n-1)
    index = 1;
    nPairs = (n-iLag)*nSamples;
    pairs = zeros(nPairs,2);
    for iSample = 1:nSamples
        for ix = 1:(n-iLag)
            pairs(index,1) = samplesMat(iSample,ix);
            pairs(index,2) = samplesMat(iSample,ix iLag);
            index = index   1;
        end
    end

    X2 = pairs(:,1)-muHat;
    Y2 = pairs(:,2)-muHat;
    ACF(iLag 1) = sum(X2.*Y2)/n*nSamples % calculate covariance
end
ACF = ACF/ACF(1); % divide by variance

Для одного временного ряда nSamples = 1 этот код выдает тот же результат, что и, autocorr несмотря на то, что он значительно медленнее. Есть ли способ использовать алгоритм FFT для вычисления ACF для нескольких временных рядов?

Комментарии:

1. Разве вы не можете просто объединить несколько реализаций и позволить autocorr выполнить усреднение внутри?

2. @AhmedFasih любая конкатенация или получение средней корреляции не дают тех же результатов, что и вычисление корреляции для матрицы пар, содержащей все пары, разделенные лагом k. Это близкое, но не то же значение. Корреляция между двумя векторами — это не то же самое, что разделение этих векторов с вычислением корреляции для каждого подразделения и усреднения.

3. Я предложил это, просто продолжая определение автокорреляции: корреляция между выборкой случайного процесса (нулевое среднее единичное отклонение) и его i обратными выборками значений E[x(t) * x(t i)] равна . Я подумал, есть ли у вас одна или десять реализаций процесса, это определение не изменится. ‘

4. Существует ряд исправлений, которые необходимо применить (см. autocorr doc и wiki ). Вам нужно вычесть среднее значение выборки внутри математического ожидания, выбрать знаменатель для нормализации суммы к среднему и т.д. Я ожидаю, что различия, которые вы обнаружите между этими методами, связаны с несоответствием этим шагам бухгалтерского учета?

5. @Greenparker к сожалению, нет бумажных ссылок. Я выяснил, как это сделать с помощью FFT, поэтому я отвечу на свой собственный вопрос.

Ответ №1:

Тот же результат получается путем вычисления среднего арифметического в области Фурье

 function [ acf ] = corrFFT(samplesMat,denomType,mu)
% Adaptation of inbuilt function autocorr to allow for multiple
% realisations, known mean and variable denominator.
%
% Input:
% samplesMat - nx by nSamples matrix of realisations
% denomType (optional) - 'cnst' or 'var'
% mu (optional) - if true mean is known

if nargin < 2
    denomType = 'cnst';
end
if nargin < 3
    mu = mean(mean(samplesMat));
end

[nx,~] = size(samplesMat);

nFFT = 2^(nextpow2(nx) 1);
F = fft(samplesMat-mu,nFFT);
F = F.*conj(F);
acf = ifft(mean(F,2));
acf = acf(1:nx);
if strcmp(denomType,'var')
    acf = (nx./[nx:-1:1]').*acf;
end

acf = acf./acf(1);
acf = real(acf);
end