1

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 },
        };
Ganessa
  • 782
  • 2
  • 7
  • 24
  • _"What did I do wrong?"_ - Without looking too deep: I _think_ you did not do anything "wrong". If the results match only by a static factor, this may well be due to normalization. Is there a difference in _that_ step that could result in a 1.5 factor difference? – Fildor Jul 17 '23 at 07:28
  • I see in the C# code you are using the fftWindowSize to normalize. So if you cahnge that, you will also change the scale of other calculations.... (not drawing anything from that, just stating the fact for now.) – Fildor Jul 17 '23 at 07:30
  • @Fildor-standswithMods Appreciate your comments. I too thought that if there was a problem, it was normalization. But excluding that process would multiply the numbers by thousands. If I exclude that process, I would have to change some more process. – Ganessa Jul 18 '23 at 00:21
  • @Fildor-standswithMods I am also concerned about that, the method of scaling to a power of 2 and filling in with zeros is described on various sites. should it be Const.NumOfSampling rather than fftWindowSize that is divided when determining the normalizedPowerSpectrum? 600Hz x 3sec = 1800 1800 < 2048 2048 / 1800 = 1.137777777777778 It differed from 1.5x. I can tell from your two comments that you completely understand what I am struggling with. Thank you. – Ganessa Jul 18 '23 at 00:26

0 Answers0