Search code examples
pythonnumpysignal-processingfft

How to extract dominant frequency from NumPy array?

I am working with simulations of populations of neurons, and would like to extract the dominant frequency from their local field potential. This takes the form of just a single vector of values, which I have plotted here. Clearly, there is some oscillatory activity happening, and I am interested in what frequency the large scale oscillations are occurring at.

Plot of LFP for a simulation

I am not very familiar with Fourier transforms outside of the general idea of them, but I believe they are the tool needed here? The sampling distance in this simulation was 0.00006 seconds or the sampling rate is 1/0.00006 Hz. I tried using the numpy.fft.fft function, but I am a little unsure about how to interpret the results. I expected a large peak around 1, corresponding to what appears like 1 Hz oscillations. Here is the code I ran, and the resulting plot.

raw = np.loadtxt(r"./file.txt")

#calculate the frequencies associated with our signal
freqs = fft.rfftfreq(len(raw),0.0006)
sp = fft.rfft(raw)

plt.plot(freqs, sp.real, freqs, sp.imag)

Entire FFT plot

I believe the high sampling rate is causing the calculation of such high frequencies, and I know I am just interested in those from ~0-4 Hz so I can plot those as:

0-4 Hz FFT plot

Or to make it easier to visualize, just ~1-4 Hz:

~1-4 Hz FFT plot

Note, I am not calculating the FFTs any differently, just plotting different regions of the results. My issue is I am unclear about how to interpret the results or also how to realize if I am making an error somewhere in my code. I can see some spikes in the results, but they are not where I expected them to be, nor am I sure on how to determine when they are significant.

Any help with this is appreciated. :^)

Solution

  • Why you don't find a peak at 1

    Just to reproduce the probable problem with my own [mre]

    import numpy as np
import matplotlib.pyplot as plt
# time with a 0.00006 second per sample
t=np.arange(0,12,0.00006)
# y is roughly sin(2π.f.t)⁴. But with some noise.
# Because of power 4, it is a 2f periodic signal. So I choose f 
# around 0.5 to mimic the 1Hz signal. And I choose not exactly 0.5
# but a more specific 0.47, to verify my ability to find it back
y=(np.sin(2*np.pi*0.47*t)+np.random.normal(0,0.1,(len(t),)))**4 

plt.plot(t,y)
plt.show()

    enter image description here

    Not the same signal, but close enough for demo purpose.

    Now, what you did is

    freqs = np.fft.rfftfreq(len(y),0.0006)
sp = np.fft.rfft(y)

    Which gives (once zoomed on the interesting part) enter image description here

    With a big value at 0, because signal is not 0 on average (it is the constant term, or 0Hz)

    Another one just before 0.1 (0.094, we surmise). And a small one (first harmonic) before 0.2.

    And then, because you were sure that anything before 1 should be ignored (but why? even if there was no error, if fft tells your the biggest peak is at 0.1 Hz, it is worth wondering why) , you zoomed on the big flat 0 area. And because of the noise, with a zoom, nothing is really flat.

    What you should have done is

    freqs = np.fft.rfftfreq(len(y),0.00006)

    Note that this has no influence whatsover on fft. The fft line doesn't even use freqs as parameter. It is only when plotting that it just changes what is indicated on x axis. Same exact plot. With just a different ticks on x-axis.

    enter image description here

    (I did not put lot of efforts to zoom exactly likewise, since I did that with the mouse. But you see that it is the same, except than what was 0.1 is now 1)

    A more accurate estimation

    You may have noted that all we see is that peak in fft is just slightly under 1. Probably my 0.94. But hard to really tell if it is 0.94, 0.96 ... Because it is a very low resolution in low frequencies. (First frequency is 0.08333... that is a period of 12 second, the the entire span. Second is twice that, so period of 6 seconds. Etc. There is no frequency in freqs between 0.91666 and 1. Hard to find accurately the 0.94 from that).

    Here, I feel that you are more looking for a repetition period than a frequency (same thing, you may say. But in fft, focus is on providing a base of periodic function that can reproduce the signal. Not to map all frequencies. For high frequencies that is almost the same thing, since resolution is High there. But for low frequencies/high periods, it is not).

    So I would rather just try to find the best autocorrelation here. For example, simply with a convolution

    conv=np.convolve(y,y[::-1])

    It gives you this enter image description here

    What is interesting, of course here, are the peaks. They occur for the shift when signals match best with itself (Think of the signal printed on 2 transparent graph. That you slide on each other).

    And there, we have not fixed a base frequency for a new vector base. Unit are "samples". What I see roughly is the position of 2nd peak, at 17760. (The one at the middle is when y is matched with itself with no shift at all. Then each peak is with a shift of 1,2,3,... or -1,-2,-3... periods) And 17760*0.06 is 1.0596 seconds. So 1/1.0596 = 0.94 Hz...