Проблема: реализация функциональности pwelch MATLAB на C#. ⇐ C#

[Anonymous]

Я пытаюсь реализовать функцию pwelch MATLAB на C#, но результаты расчета PSD отличаются. Функция MATLAB используется следующим образом:
[fzz, freq] = pwelch(data, hanning(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 windowList = new List();
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;
}
}

Почему я не получаю таких же результатов?

Подробнее здесь: https://stackoverflow.com/questions/792 ... in-c-sharp