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 },
};