#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