Matlab中使用FFT的多个时间序列的自相关答案

【问题标题】：Autocorrelation of multiple time series in Matlab using FFTMatlab中使用FFT的多个时间序列的自相关
【发布时间】：2017-02-21 08:46:14
【问题描述】：

Matlab 函数autocorr（文档here）使用快速傅里叶变换 (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 by 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？

【问题讨论】：

你不能把多个实现串联起来，让autocorr在里面做平均吗？
@AhmedFasih 任何类型的连接或采用平均相关性都不会给出与计算包含所有由滞后 k 分隔的对的对矩阵的相关性相同的结果。它很接近，但值不同。两个向量之间的相关性与对这些向量进行细分计算每个细分的相关性并取平均值不同。
我建议这只是关于自相关的定义：（零均值单位方差）随机过程的样本与其值i 样本之间的相关性是E[x(t) * x(t+i)]。我在想，无论你有一个或十个过程的实现，这个定义都不会改变。'
需要进行一些更正（参见autocorr doc 和wiki）。您需要在期望值中减去样本均值，选择用于将总和归一化为均值的分母等。我希望您发现这些方法之间的差异来自不匹配这些簿记步骤？
@Greenparker 没有论文参考对不起。我确实想出了如何用 FFT 来做，所以我会回答我自己的问题。

标签： matlab time-series fft correlation

【解决方案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

【讨论】：