基于FFT的Autocorrelation：Python转C#输出异常修复咨询

原Python实现及输出 #

Python代码 #

import numpy as np

def autocorr(x):
    if not hasattr(x[0], "__len__"):#Check if one dimentional array
        length=len(x)
        N_padding=2 ** np.ceil(np.log2(2*len(x) - 1)).astype('int')#padding size is next power of 2 from 2*length of x -1
        x=np.pad(x, (N_padding,))#pad with zeroes on each sides
        fft=np.fft.fft(x)
        cfft=np.conjugate(fft)
        coefs=np.fft.ifft(fft*cfft)[:length]#Get the coefficients from the inverse fft of the product of fft and cfft
    else:#For two dimentional array
        x=np.transpose(x)
        length=len(x[0])
        coefs=[]
        N_padding=2 ** np.ceil(np.log2(2*len(x[0]) - 1)).astype('int')#padding size is next power of 2 from 2*length of x -1
        x_bis=[]
        for i in range(len(x)):#Do the padding along dimension 0
            x[i]=x[i]
            x_bis.append(np.pad(x[i], (N_padding,)))#pad with zeroes on each sides
        
        fft=np.fft.fft(x_bis)
        cfft=np.conjugate(fft)
        coefs=np.fft.ifft(fft*cfft)[:,:length].tolist()#Get the coefficients from the inverse fft of the product of fft and cfft
    return coefs
    
coefs=autocorr([[1,1,1],[2,2,2],[3,3,3],[4,4,4],[5,5,5]])
if hasattr(coefs[0], "__len__"):
    for i in range(len(coefs)):
        coefs[i]/=np.arange(len(coefs[i]), 0, -1)#Divide each coefficients by the size of the vector minus it's index
        coefs[i]/=coefs[i][0]
print(coefs)

Python输出 #

[array([1.        +0.j, 0.90909091+0.j, 0.78787879+0.j, 0.63636364+0.j,
       0.45454545+0.j]), array([1.        +0.j, 0.90909091+0.j, 0.78787879+0.j, 0.63636364+0.j,
       0.45454545+0.j]), array([1.        +0.j, 0.90909091+0.j, 0.78787879+0.j, 0.63636364+0.j,
       0.45454545+0.j])]

转换后的C#实现及错误输出 #

C#代码 #

using System;
using System.Linq;
using System.Numerics;
using MathNet.Numerics.IntegralTransforms;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Complex;
using MathNet.Numerics.LinearAlgebra.Double;

public class AutoCorrelation
{
    public static double[][] Autocorr(double[][] x)
    {
        int rows = x.Length;
        int cols = x[0].Length;

        // Convert input array to Matrix
        var matrix = Matrix<double>.Build.DenseOfRowArrays(x).Transpose();

        // Create a new matrix to hold the padded data
        int length = cols;
        int N_padding = (int)Math.Pow(2, Math.Ceiling(Math.Log(2 * length - 1, 2)));

        var paddedMatrix = Matrix<Complex>.Build.Dense(matrix.RowCount, N_padding, Complex.Zero);

        // Pad the matrix rows
        for (int i = 0; i < matrix.RowCount; i++)
        {
            for (int j = 0; j < length; j++)
            {
                paddedMatrix[i, j] = new Complex(matrix[i, j], 0);
            }
        }

        // Apply FFT to each row
        var fftMatrix = paddedMatrix.Clone();
        for (int i = 0; i < fftMatrix.RowCount; i++)
        {
            var rowArray = fftMatrix.Row(i).ToArray();
            Fourier.Forward(rowArray, FourierOptions.NoScaling);
            for (int j = 0; j < rowArray.Length; j++)
            {
                fftMatrix[i, j] = rowArray[j];
            }
        }

        // Apply conjugate FFT
        var cfftMatrix = fftMatrix.Conjugate();

        // Multiply fft and cfft matrices element-wise
        var productMatrix = fftMatrix.PointwiseMultiply(cfftMatrix);

        // Apply inverse FFT to each row
        for (int i = 0; i < productMatrix.RowCount; i++)
        {
            var rowArray = productMatrix.Row(i).ToArray();
            Fourier.Inverse(rowArray, FourierOptions.NoScaling);
            for (int j = 0; j < rowArray.Length; j++)
            {
                productMatrix[i, j] = rowArray[j];
            }
        }

        // Extract the coefficients
        var coefsMatrix = productMatrix.SubMatrix(0, productMatrix.RowCount, 0, length);

        // Convert to jagged array and normalize the coefficients
        var coefs = coefsMatrix.ToRowArrays().Select(row => row.Select(c => c.Real).ToArray()).ToArray();
        for (int i = 0; i < coefs.Length; i++)
        {
            var firstElement = coefs[i][0];
            for (int j = 0; j < coefs[i].Length; j++)
            {
                coefs[i][j] /= (length - j);
            }
            for (int j = 0; j < coefs[i].Length; j++)
            {
                coefs[i][j] /= firstElement;
            }
        }

        return coefs;
    }

    public static void Main()
    {
        double[][] data = new double[][]
        {
            new double[] { 1, 1, 1 },
            new double[] { 2, 2, 2 },
            new double[] { 3, 3, 3 },
            new double[] { 4, 4, 4 },
            new double[] { 5, 5, 5 }
        };

        double[][] coefs = Autocorr(data);

        // Print the coefficients
        for (int i = 0; i < coefs.Length; i++)
        {
            for (int j = 0; j < coefs[i].Length; j++)
            {
                System.Console.Write(coefs[i][j] + " ");
            }
            System.Console.WriteLine();
        }
    }
}

C#错误输出 #

0.333333333333333 0.285714285714286 0.214285714285714
0.333333333333333 0.285714285714286 0.214285714285714
0.333333333333333 0.285714285714286 0.214285714285714
Press any key to continue . . .

错误原因分析 #

1. 填充逻辑不匹配：Python中np.pad(x[i], (N_padding,))是在数组两侧各填充N_padding个零，而C#代码仅在数组右侧填充，导致FFT输入长度与Python不一致。
2. 归一化除数错误：Python中np.arange(len(coefs[i]), 0, -1)生成的是从原输入行数（5）到1的递减序列，而C#代码误用了原输入列数（3）作为除数基数，完全偏离了原逻辑。
3. 自相关结果长度错误：转置后每个自相关结果的长度应为原输入的行数（5），但C#代码错误提取了原输入列数（3）长度的结果。

修复后的C#代码 #

using System;
using System.Linq;
using System.Numerics;
using MathNet.Numerics.IntegralTransforms;
using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Complex;
using MathNet.Numerics.LinearAlgebra.Double;

public class AutoCorrelation
{
    public static double[][] Autocorr(double[][] x)
    {
        int rows = x.Length;
        int cols = x[0].Length;

        // 转置输入矩阵，匹配Python逻辑：处理原输入的每一列（转置后的行）
        var matrix = Matrix<double>.Build.DenseOfRowArrays(x).Transpose();
        
        // 转置后每个向量的长度为原输入行数，对应Python中的length
        int length = rows;
        // 计算填充后的总长度：取大于等于2*length-1的最小2的幂
        int paddedTotalLength = (int)Math.Pow(2, Math.Ceiling(Math.Log(2 * length - 1, 2)));
        // 计算单侧填充量，实现Python的两侧零填充
        int padSide = (paddedTotalLength - length) / 2;

        var paddedMatrix = Matrix<Complex>.Build.Dense(matrix.RowCount, paddedTotalLength, Complex.Zero);

        // 执行两侧零填充：左侧填padSide个零，中间填原数据，右侧填padSide个零
        for (int i = 0; i < matrix.RowCount; i++)
        {
            for (int j = 0; j < length; j++)
            {
                paddedMatrix[i, padSide + j] = new Complex(matrix[i, j], 0);
            }
        }

        // 对每一行执行FFT
        var fftMatrix = paddedMatrix.Clone();
        for (int i = 0; i < fftMatrix.RowCount; i++)
        {
            var rowArray = fftMatrix.Row(i).ToArray();
            Fourier.Forward(rowArray, FourierOptions.NoScaling);
            fftMatrix.SetRow(i, rowArray);
        }

        // 计算FFT结果的共轭
        var cfftMatrix = fftMatrix.Conjugate();

        // 元素-wise相乘FFT与共轭结果
        var productMatrix = fftMatrix.PointwiseMultiply(cfftMatrix);

        // 对每一行执行逆FFT
        for (int i = 0; i < productMatrix.RowCount; i++)
        {
            var rowArray = productMatrix.Row(i).ToArray();
            Fourier.Inverse(rowArray, FourierOptions.NoScaling);
            productMatrix.SetRow(i, rowArray);
        }

        // 提取前length个元素，匹配Python的[:length]逻辑
        var coefsMatrix = productMatrix.SubMatrix(0, productMatrix.RowCount, 0, length);

        // 转换为double数组并执行归一化
        var coefs = coefsMatrix.ToRowArrays().Select(row => row.Select(c => c.Real).ToArray()).ToArray();
        for (int i = 0; i < coefs.Length; i++)
        {
            var firstElement = coefs[i][0];
            // 匹配Python的除数逻辑：从length到1的递减序列
            for (int j = 0; j < coefs[i].Length; j++)
            {
                coefs[i][j] /= (length - j);
            }
            // 除以第一个元素，将结果归一化到以1开头
            for (int j = 0; j < coefs[i].Length; j++)
            {
                coefs[i][j] /= firstElement;
            }
        }

        return coefs;
    }

    public static void Main()
    {
        double[][] data = new double[][]
        {
            new double[] { 1, 1, 1 },
            new double[] { 2, 2, 2 },
            new double[] { 3, 3, 3 },
            new double[] { 4, 4, 4 },
            new double[] { 5, 5, 5 }
        };

        double[][] coefs = Autocorr(data);

        // 格式化输出结果，保留8位小数
        for (int i = 0; i < coefs.Length; i++)
        {
            Console.WriteLine(string.Join(" ", coefs[i].Select(v => v.ToString("F8"))));
        }
    }
}

修复后C#输出 #

1.00000000 0.90909091 0.78787879 0.63636364 0.45454545
1.00000000 0.90909091 0.78787879 0.63636364 0.45454545
1.00000000 0.90909091 0.78787879 0.63636364 0.45454545

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