Reproduce the results of using Python's spectrogram() in C#

54 Views Asked by At

I would like you to point out various insights from the perspective of a C# expert, a Python expert, and an EEG expert.

I would like to get the same results using MathNet in C# as I get using scipy.signal.spectrogram in Python.

https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.spectrogram.html https://numerics.mathdotnet.com/api/MathNet.Numerics.IntegralTransforms/Fourier.htm

I got to the point where I could get pretty close data. I multiply the C# results by 1.5 to get results that are about equal to the Python results. But I do not want to use a magic number. I want to find out why and fix my code. Could you please point out the errors in my code?

The fact that they match by a factor of 1.5 leads me to believe I am making a mistake in multiplication or division. I moved the timing of using the correction factor for the Hamming window back and forth, which produced more different results.

What did I do wrong?

Python

# Sampling rate [Hz]
fs = 600.0

# FFT window length [sec].
win_sec = 3.0

# Upper frequency [Hz]
fq_lim = 50

# Sampling interval [sec]
dt = 1. / fs

# Data length [points]
n = int(fs * win_sec)

# Power of 2 greater than argument calculated
nfft = 2 ** math.ceil(math.log2(n))

# Frequency signal interval [/sec].
df_f = 1.0 / (nfft * dt)

# Spectrogram calculation Frequency (freq) Time (t) Power spectral density PSD (μV²/Hz) calculation
freq, t, spctgram = signal.spectrogram(sig, fs=fs, nperseg=n, nfft=nfft, noverlap=(win_sec-1)*fs, window='hamming', mode='psd', return_onesided=True)

C#

        public Dictionary<string, float> Analyze(List<float> brainwaveData)
        {
            // power of 2 greater than the sampling number
            var fftWindowSize = (int)Math.Pow(2, Math.Ceiling(Math.Log(Const.NumOfSampling, 2)));

            // Fill with zeros until sampling number is next power of 2
            while (brainwaveData.Count < fftWindowSize)
            {
                brainwaveData.Add(0f);
            }

            // Apply Hamming window for FFT
            var coefficientHammingWindow = 0.54f;
            double[] window = Window.HammingPeriodic(fftWindowSize);
            var signalWindow = brainwaveData.Select((d, i) => new Complex32((float)(d * window[i] / coefficientHammingWindow), 0)).ToArray();

            // Execute FFT
            Fourier.Forward(signalWindow, FourierOptions.Matlab);

            // Power spectrum [μV²] calculation Amplitude squared
            var powerSpectrum = signalWindow.Select(d => (d.Real * d.Real) + (d.Imaginary * d.Imaginary)).ToArray();

            // power spectrum normalization [μV²]
            var normalizedPowerSpectrum = powerSpectrum.Select(d => d / fftWindowSize).ToArray();

            // Power spectral density [μV²/Hz] calculation
            var powerSpectrumDensity = normalizedPowerSpectrum.Select(d => d / Const.SamplingFrequencyHz).ToArray();

            // Calculated for each EEG component
            Dictionary<string, float> result = new Dictionary<string, float>();
            foreach (var listBand in listBands)
            {
                result.Add(listBand.BandName, calcFrequencyBandAverage(powerSpectrumDensity, Const.SamplingFrequencyHz, listBand.MinFreqHz, listBand.MaxFreqHz));
            }
            return result;
        }

        private float calcFrequencyBandAverage(float[] spectrum, int sampleRate, double minHz, double maxHz)
        {
            // Find the index for a given frequency band (Hz)
            var minIndex = (int)Math.Ceiling(spectrum.Length * minHz / sampleRate);
            var maxIndex = (int)Math.Floor(spectrum.Length * maxHz / sampleRate);

            // Calculate average
            return spectrum.Skip(minIndex).Take(maxIndex - minIndex).Average();
        }

        List<BrainwaveFrequencyBands> listBands = new List<BrainwaveFrequencyBands>
        {
            new BrainwaveFrequencyBands{ BandName = "Delta", MinFreqHz = 1f, MaxFreqHz = 4f },
            new BrainwaveFrequencyBands{ BandName = "Theta", MinFreqHz = 4f, MaxFreqHz = 8f },
            new BrainwaveFrequencyBands{ BandName = "Alpha1", MinFreqHz = 8f, MaxFreqHz = 10f },
            new BrainwaveFrequencyBands{ BandName = "Alpha2", MinFreqHz = 10f, MaxFreqHz = 12f },
            new BrainwaveFrequencyBands{ BandName = "Beta1", MinFreqHz = 12f, MaxFreqHz = 20f },
            new BrainwaveFrequencyBands{ BandName = "Beta2", MinFreqHz = 20f, MaxFreqHz = 30f },
            new BrainwaveFrequencyBands{ BandName = "Gamma1", MinFreqHz = 30f, MaxFreqHz = 40f },
            new BrainwaveFrequencyBands{ BandName = "Gamma2", MinFreqHz = 40f, MaxFreqHz = 50f },
        };
0

There are 0 best solutions below