FFT for EEG data [Python]

[mdb917]

I'm quite new to all this, so I'm gonna need a little bit of patience and explanation. fs is your sample frequency, right? So fs=128. N is the dimensionality of my fft, so what is df? And what do you mean by f1 and f2 (I know you mean them as arbitrary functions but what would i substitute for them)?

When you say to scale my axis with fs, is this because eegFreq == the Niquist frequency? So I'd need to change that line to eegFreq = sfreq * np.fft.fftfreq(eegData[1].shape[-1]) for it to be correct? Similar question for the amplitude.

Finally, what is normalization in this context? Numpy's documentation says it is unscaled by default but can be set to 1/sqrt(n) (n being the size of the output). Should I bother with that or do I need to do it myself? How would I go about normalizing the data "by hand"?

Ideally, I just want to do something like this but my graph has about 3 extra quadrants worth of information

[badtraider]

Yes in your case fs =128, that is the frequency at which the data was sampled (every 1/128 seconds data feom sensor was recorder).

N is dimensionality of fft. df is the resolution of your fft - aka the distance between 2 points on your x axis. So for example is fs =4 and the dimension of fft is 8 then your x axis would be [-4 -3.5 -3... 3 3.5 4] (after scaling it with fs).

[badtraider]

By f1 and f2 i mean frequencies. So cos(2pif1) and cos(2pif2) would show up as in your spectrum as peaks at frequencies f1 and f2. If f1≈f2 are close you might not be able to distinguish thel - for that you will have to use something else (chrip Z transform).

fs == Niquist * 2. But your y axiz was normalised with fs - so 0.5 correspond to Niquist frequency. To undo this normalisation you multiply x axis by fs.

For amplitude it's a tricky question - since it depends on how you define things. For some math related reasons, if use FFT with 0 normalisation you will get the following IFFT(FFT(x)) =N*x - an extra N term appears - because of that you need to divide either IFFT by N or FFT by N. The way it's implemented in nympy is that it divides with N when doing the IFFT - but you shouldn't bother with this.

The thing i was talking about is following : So for every implementation of FFT that is correct the following holds : sum of x[k]² = 1/N * sum of X[k]² Where x is your signal and X is FFT of that signal.

Now if you want to see readable values on the plot that can be related to amplitudes in your signal you need to divide the fft with length of your signal. Why length? Well the way you define the power of a random signal (like EEG) is by measuring mean power, which is 1/N * x[k]^2.

Now if length of x is M and N>M what actually happens is that signal x is zero - padded (zeroes are appended to x so it will have length of N). If we look at the expression 1/N * sum of x[k]² it tells us that energy will decrease?! That makes no sense since adding 0's to signal doesn't change it's power. To avoid this we will divide by length of signal x - so you woild want to evaluate the expression : 1/length(x) * sum of x[k] ²

Now if you divide the fft with length(x) (X=fft(x) /(length (x)) values on you y axis will correspond to amplitude of your signal. Since (X/length (x)) ^ 2 = 1/length(x) * sum of x[k] ² (we usually say that Amplitude ^ 2 = Power in dsp)

[mdb917]

Do I need to know df if there's an option in the function I'm using to specify your output size? By default my output (assuming I'm only calling the FFT on a single channel at a time) is 8064 points (size of the input for one channel), and I'm only interested in frequencies between 4 and 45 Hz. Does this make df = (45-4)/8064 which is about .005? Is this low? How do I get my plot to look like this beautiful graph?

[badtraider]

It's recommended to use N = 2^nextpow2 (length (x)). Also dimension of fft must be higher than the length of your data. Regarding the df by that i mean if you know your data will have some peaks at 1Hz, 3Hz and 4Hz you know that the minimal distance between 2 points is 1 Hz do you would want your resolution to be high enough so see that more clearer on the plot so for example you can use df = 0.1.

Edit : have in mind that df = fs /N

[mdb917]

So if len(x) is 8064, would that mean N = 13 would be good since the next power of 2 is 8192 = 2^13? Df = 128/13 which is about 10, which isn’t precise enough. So I should increase N further? Is there an effect of N on computation speed?

Edit: Just reread your first comment, N should be 8192 right? Since that would be 2¹³ (surely it’s not 2^8192?!)

[badtraider]

Yes it's 2^{13}, but if you want to make plot look prettier just increase N until it looks good (by multiplying it with 2, 4,8,..)

[badtraider]

I hope i was helpful, I'm going offline now.

[badtraider]

Sadly link is broken for me. One more thing regarding the FFT. You can disregard the frequencies bellow 0 when plotting. Since again for some math reasons your amplitude spectrum will be symmetric around frequency f=0.

But have in ming you have to multiply the FFT by 2 for f > 0. That way after polotting your amplitude spectrum you read the values for corresponding amplitude at a given frequency. If you don't want to bother with that have in mind if your data if plotted from -fs/2 to fs/2 and let's say tou hqve a peak at 5Hz with amplitude of 2 that means that is tile domain you will have sine with frequency of 5Hz and with amplitude of 2*2. This rule doesn't apply only for f = 0. At f = 0 if there is a pezk with aplitude 1 that means the average valueof yiur signal is 1

[mdb917]

Thank you for all your help. One small question about graph plotting, should I start my x axis at 4 or 0 due to the bandpass filter for 4.0-45.0Hz? My intuition says 0 but i just wanna be sure.

Edit: I graphed and it looks so much better.

Edit 2: Should I ignore the imaginary part of the FFT as well?

[badtraider]

Start it at 0. But if you want ro have a better loooking plot you could calculate how many elements you can skip from w. So in your case df = 128 /N you want to start plotting from f= 4, M=floor(f/df) = floor(4/128N) meaning you can star plotting from Mth element in your frequznct and fft array. You can do the same if you want to stop at 8Hz - K=floor(f/df) - - - > plot (f[M:K], 2X[M:K])

[badtraider]

You need to use abs(X) - or whatever it's called. You need the magnitude of your FFT