如何修正基于MathNet.Numerics的3D向量FFT自相关计算代码？

核心问题分析 #

原FFT实现的错误在于将3D向量转换为模值后计算标量自相关，这和原循环中逐对向量点积求和再平均的逻辑完全不符。正确的数学逻辑是：3D向量序列的自相关（延迟τ处的值）等于三个分量各自自相关（延迟τ处的值）的总和，再除以有效样本数。

修正步骤 #

1. 将3D向量序列拆解为X、Y、Z三个独立的标量序列
2. 对每个标量序列分别用FFT计算自相关
3. 将三个分量的自相关结果按延迟τ对应相加，得到点积的自相关总和
4. 对每个延迟τ，用总和除以有效样本数（总样本数 - τ）得到最终平均值

代码实现（C# + MathNet.Numerics） #

using MathNet.Numerics;
using MathNet.Numerics.IntegralTransforms;
using System.Numerics;
using System.Collections.Generic;

public class Vector3AutoCorrelation
{
    public double[] ComputeAutoCorrelationFFT(List<Vector3> vectors)
    {
        int n = vectors.Count;
        if (n == 0) return Array.Empty<double>();

        // 拆解三个分量序列
        double[] x = vectors.Select(v => v.X).ToArray();
        double[] y = vectors.Select(v => v.Y).ToArray();
        double[] z = vectors.Select(v => v.Z).ToArray();

        // 计算每个分量的FFT自相关
        double[] xCorr = ComputeScalarAutoCorrelationFFT(x);
        double[] yCorr = ComputeScalarAutoCorrelationFFT(y);
        double[] zCorr = ComputeScalarAutoCorrelationFFT(z);

        // 合并三个分量的自相关，得到点积的自相关总和并归一化
        double[] result = new double[n];
        for (int τ = 0; τ < n; τ++)
        {
            result[τ] = (xCorr[τ] + yCorr[τ] + zCorr[τ]) / (n - τ);
        }

        return result;
    }

    private double[] ComputeScalarAutoCorrelationFFT(double[] sequence)
    {
        int n = sequence.Length;
        int fftSize = NextPowerOfTwo(n * 2 - 1); // 避免循环卷积混叠

        // 补零到FFT尺寸
        Complex[] input = new Complex[fftSize];
        for (int i = 0; i < n; i++)
        {
            input[i] = new Complex(sequence[i], 0);
        }

        // FFT正变换
        Fourier.Forward(input, FourierOptions.Matlab);

        // 计算功率谱（共轭相乘）
        Complex[] powerSpectrum = input.Select(c => c * Complex.Conjugate(c)).ToArray();

        // 逆FFT得到自相关（未归一化）
        Fourier.Inverse(powerSpectrum, FourierOptions.Matlab);

        // 提取前n个点，取实部（虚部为浮点误差可忽略）
        double[] correlation = new double[n];
        for (int τ = 0; τ < n; τ++)
        {
            correlation[τ] = powerSpectrum[τ].Real;
        }

        return correlation;
    }

    private int NextPowerOfTwo(int value)
    {
        int power = 1;
        while (power < value) power <<= 1;
        return power;
    }
}

关键说明 #

• 补零到2n-1的下一个2的幂，是为了避免循环卷积的混叠问题，确保线性卷积结果正确
• 最终除以n-τ是对齐原循环逻辑：延迟τ时仅存在n-τ个有效点积对，保证平均值计算一致
• 该实现时间复杂度为O(n log n)，对比原循环的O(n²)，处理千万级数据时性能会有数量级提升

正确性验证 #

用小规模数据集同时运行原循环实现和FFT实现，对比各延迟τ的结果，误差应在浮点精度范围内（可忽略）。

内容的提问来源于stack exchange，提问作者user366312