如何在矢量数据的情况下应用FFT？

Question

我编写了以下类来计算 3d 向量列表的自相关。

我从这个链接

获取了公式

public static class AutocorrVec3
{
    private static double C(int t, List<Vec3> vectors)
    {
        int n = vectors.Count;
        if (t >= n || t < 0)
            throw new ArgumentException("Invalid value for t. It must be between 0 and n-1.");

        double sum = 0;
        for (int i = 0; i < n - t; i++)
        {
            sum += vectors[i].Dot(vectors[i + t]);
        }

        return sum / (n - t);
    }

    public static (List<double> taus, List<double> autocorrs) 
        GetAutoCorrelationPoints(List<Vec3> vectors, int maxLag)
    {
        var tValues = new List<double>();
        var cResults = new List<double>();

        double c0 = C(0, vectors); // This is the normalization factor
        Console.WriteLine($"Normalization factor: {c0}");

        for (int lag = 0; lag <= maxLag; lag++)  // Start from 0 to include the autocorrelation at lag 0
        {
            try
            {
                double cValue = C(lag, vectors);

                Console.WriteLine($"Lag={lag}, Raw Autocorr={cValue}, Normalized Autocorr={cValue / c0}");
                tValues.Add(lag);
                cResults.Add(cValue / c0); // Normalization is done here
            }
            catch (ArgumentException ex)
            {
                Console.WriteLine(ex.Message);
                break;
            }
        }
        return (tValues, cResults);
    }
}

问题是，

GetAutoCorrelationPoints()

非常慢。例如，我需要转换 24 个文件，每个文件包含 10000000 个 3d 矢量。 24小时后，连一个数据文件都无法完成计算。

在这种情况下我该如何应用FFT？

我想用

MathNet.Numerics

。

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


namespace FourierTransformCSharp
{
    public static class AutocorrVec3
    {
        // Compute the autocorrelation of a time series using FFT
        public static double[] ComputeAutocorrelationUsingFFT(List<Vec3> vectors)
        {
            int n = vectors.Count;
            // Create a zero-padded list double the size of the original for FFT
            var paddedVectors = new Complex[2 * n];
            for (int i = 0; i < n; i++)
            {
                // Convert vector to complex number with magnitude as real part
                paddedVectors[i] = new Complex(vectors[i].Magnitude(), 0);
            }
            for (int i = n; i < 2 * n; i++) // Zero padding
            {
                paddedVectors[i] = Complex.Zero;
            }

            // Perform FFT on the zero-padded list
            Fourier.Forward(paddedVectors, FourierOptions.Default);

            // Compute power spectrum (magnitude squared)
            for (int i = 0; i < paddedVectors.Length; i++)
            {
                var magnitude = paddedVectors[i].Magnitude;
                paddedVectors[i] = new Complex(magnitude * magnitude, 0);
            }

            // Perform Inverse FFT to obtain the autocorrelation function
            Fourier.Inverse(paddedVectors, FourierOptions.Default);

            // Extract the real parts as the autocorrelation result
            var autocorr = new double[n];
            for (int i = 0; i < n; i++)
            {
                autocorr[i] = paddedVectors[i].Real;
            }

            // Normalize the autocorrelation result
            var normalizationFactor = autocorr[0];
            for (int i = 0; i < n; i++)
            {
                autocorr[i] /= normalizationFactor;
            }

            return autocorr;
        }

        // Calculate autocorrelation of vector time series at lag t
        public static double C(int t, List<Vec3> vectors)
        {
            double[] autocorr = ComputeAutocorrelationUsingFFT(vectors);
            return autocorr[t];
        }

        // Get autocorrelation values for lags from 0 to maxLag
        public static (List<int> taus, List<double> autocorrs) GetAutoCorrelationPoints(List<Vec3> vectors, int maxLag)
        {
            List<int> taus = new List<int>();
            List<double> autocorrs = new List<double>();
            double[] autocorr = ComputeAutocorrelationUsingFFT(vectors);

            for (int t = 0; t <= maxLag && t < vectors.Count; t++)
            {
                taus.Add(t);
                autocorrs.Add(autocorr[t]);
            }

            return (taus, autocorrs);
        }

        // Parallel computation is unnecessary as FFT-based autocorrelation is already fast
        // Use GetAutoCorrelationPoints for parallel computation
    }

    public class Vec3
    {
        public double X { get; }
        public double Y { get; }
        public double Z { get; }

        public Vec3(double x, double y, double z)
        {
            X = x;
            Y = y;
            Z = z;
        }

        // Compute the magnitude of the vector
        public double Magnitude()
        {
            return Math.Sqrt(X * X + Y * Y + Z * Z);
        }
    }

    public class AutocorrelationDriver
    {
        public static void Main()
        {
            // Define your list of 3D vectors
            List<Vec3> vectors = new List<Vec3>
            {
                new Vec3(1, 1, 1),
                new Vec3(2, 2, 2),
                new Vec3(3, 3, 3),
                new Vec3(4, 4, 4),
                new Vec3(5, 5, 5)
            };

            // Compute the autocorrelation
            var aurocorr = AutocorrVec3.GetAutoCorrelationPoints(vectors, vectors.Count - 1);

            // Output the results
            Console.WriteLine("Lag\tAutocorrelation");
            for (int i = 0; i < aurocorr.taus.Count; i++)
            {
                Console.WriteLine($"{aurocorr.taus[i]}\t{aurocorr.autocorrs[i]}");
            }

            Console.ReadKey();
        }
    }
}

上面的代码是我写的。

正确的输出如下：

但是，我编写的代码给出了以下输出：

Lag     Autocorrelation
0       1
1       0.727272727272727
2       0.472727272727273
3       0.254545454545455
4       0.0909090909090911

如何修复我的代码？

Answer 1

就目前而言，基于 FFT 的计算实际上没有任何问题。您记得在所有正确的位置进行零填充（做得很好），它给出了我对您的测试数据所期望的答案。

这两个代码之间的差异是因为您忘记了通过贡献滞后分量 1/(n-t) 的数量来标准化 FFT 计算的输出。

应用该校正进行 FFT 相关性计算可得到与玩具测试数据上的普通方法完全相同的答案。

我不太确定将相关函数应用于矢量大小的优点，但这完全是另一回事。这是示例数据的表格。

数据	滞后	sum(a(i).a(i+t))	规范	/(n-t)	范数 C(t,v5)
1	0	55	1	0.2	1
2	1	40	0.727272727	0.181818182	0.909090909
3	2	26	0.472727273	0.157575758	0.787878788
4	3	14	0.254545455	0.127272727	0.636363636
5	4	5	0.090909091	0.090909091	0.454545455

您可以通过一种方法来进一步优化它，方法是返回向量中的幂，而不是返回magnitude = sqrt(X.X+Y.Y+Z.Z)，然后在返回“X.X+Y.Y+”时对其进行平方Z.Z" 它已准备好用于计算功率谱。我真的不确定这个计算有什么物理解释。

顺便说一句，通过使用 FFT 的实数到复数共轭对称版本，您几乎可以将速度提高一倍并将内存需求减半。这避免了将真实数据提升为零虚部的复杂数据。首先尝试“按原样”，因为对于较大的 N，Nlog N 比 N^2 小很多。

但是速度提高两倍可能仍然值得付出额外的努力。