问题：在 C# 中实现 MATLAB 的 pwelch 功能

Question

我尝试用C#实现MATLAB的pwelch函数，但PSD计算结果不同。 MATLAB函数的使用方法如下： [fzz, 频率] = pwelch(数据, 汉宁(512), 256, 512, fs);

在C#中，我实现了如下所示的功能。当我运行 [fzz, freq] = pwelch(data, hanning(1024), 512, 1024, fs); 时，结果与我的 C# 实现中的结果相同。但是，当我运行 [fzz, freq] = pwelch(data, hanning(512), 256, 512, fs); 时，结果与我的 C# 实现不同。

为了确保一致性，我从 MATLAB 中提取了窗口数据并直接在 C# 实现中使用它。

造成这种差异的原因是什么？

[C# 代码]

public class PWELCH
{
    public static (double[] frequencies, double[] psd) Welch(double[] data, double samplingRate, int windowLength, int overlap, int nfft, string fpath_windows)
    {
        // Step size
        int step = windowLength - overlap;
        int numSegments = (int)Math.Floor((double)(data.Length - overlap) / step);

        // Hanning window
        //double[] window = GenerateHanningWindow(windowLength);
        double[] window = GenerateHanningWindow(fpath_windows);

        // Frequency vector
        double[] frequencies = new double[nfft / 2 + 1];
        for (int i = 0; i < frequencies.Length; i++)
        {
            frequencies[i] = i * (samplingRate / nfft);
        }

        // PSD array
        double[] psd = new double[nfft / 2 + 1];
        double windowEnergy = 0;

        //Calculate window energy
        for (int i = 0; i < window.Length; i++)
        {
            windowEnergy += window[i] * window[i];
        }


        // Welch algorithm
        for (int i = 0; i < numSegments; i++)
        {
            // Extract segment
            int start = i * step;
            double[] segment = new double[windowLength];
            Array.Copy(data, start, segment, 0, windowLength);

            // Apply window
            for (int j = 0; j < windowLength; j++)
            {
                segment[j] *= window[j];
            }

            // Zero-pad to nfft and calculate FFT
            Complex[] fftResult = new Complex[nfft];
            for (int j = 0; j < nfft; j++)
            {
                if (j < windowLength)
                {
                    fftResult[j] = new Complex(segment[j], 0); // double 값 그대로 사용
                }
                else
                {
                    fftResult[j] = Complex.Zero;
                }
            }

            Fourier.Forward(fftResult, FourierOptions.Matlab);

            // Calculate PSD
            for (int j = 0; j < psd.Length; j++)
            {
                psd[j] += fftResult[j].MagnitudeSquared() / (windowEnergy * samplingRate);
            }

        }

        // Average over segments
        for (int i = 0; i < psd.Length; i++)
        {
            psd[i] /= numSegments;
        }

        // Apply one-sided spectrum scaling
        psd[0] /= 2; // DC component (0 Hz)
        psd[psd.Length - 1] /= 2; // Nyquist frequency

        for (int i = 1; i < psd.Length - 1; i++)
        {
            psd[i] *= 2;
        }

        return (frequencies, psd);

    }

    public static double[] GenerateHanningWindow(string fpath)
    {
        // Read all lines from the file
        string[] lines = File.ReadAllLines(fpath);

        // Convert lines to a double array
        List<double> windowList = new List<double>();
        foreach (var line in lines)
        {
            if (double.TryParse(line, out double value))
            {
                windowList.Add(value);
            }
            else
            {
                throw new FormatException($"Invalid data in file: '{line}' is not a valid double.");
            }
        }

        // Convert the list to an array and return
        return windowList.ToArray();
    }

    public static double AdjustSamplingRateToPowerOfTwo(double samplingRate)
    {
        // Calculate the nearest power of 2
        double log2Value = Math.Log(samplingRate, 2);
        int lowerPower = (int)Math.Floor(log2Value);
        int upperPower = (int)Math.Ceiling(log2Value);

        // Compare the distances to the nearest powers of 2
        double lowerValue = Math.Pow(2, lowerPower);
        double upperValue = Math.Pow(2, upperPower);

        return (samplingRate - lowerValue < upperValue - samplingRate) ? lowerValue : upperValue;
    }
}

为什么我没有得到相同的结果？

Answer 1

使用较小的窗口大小（如 hanning(512)）时，Welch 方法的 C# 实现与 MATLAB 的 pwelch 函数之间的差异可能是由于分段数据时在信号边缘处理加窗方式的差异所致；具体来说，当计算窗口大小小于段长度的段时，MATLAB 如何在边界处隐式地对信号进行零填充。