Power spectral density#

Hide code cell source
import pynapple as nap
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
custom_params = {"axes.spines.right": False, "axes.spines.top": False}
sns.set_theme(style="ticks", palette="colorblind", font_scale=1.5, rc=custom_params)

Generating a signal#

Let’s generate a dummy signal with 2Hz and 10Hz sinusoide with white noise.

F = [2, 10]

Fs = 2000
t = np.arange(0, 200, 1/Fs)
sig = nap.Tsd(
    t=t,
    d=np.cos(t*2*np.pi*F[0])+np.cos(t*2*np.pi*F[1])+2*np.random.normal(0, 3, len(t)),
    time_support = nap.IntervalSet(0, 200)
    )
Hide code cell source
plt.figure()
plt.plot(sig.get(0, 0.4))
plt.title("Signal")
plt.xlabel("Time (s)")
Text(0.5, 0, 'Time (s)')
../_images/e0289ffcb937c602bdbc862ff78053ff802eae823833eebd6cec5029842f90b5.png

Computing FFT#

To compute the FFT of a signal, you can use the function nap.compute_fft. With norm=True, the output of the FFT is divided by the length of the signal.

fft = nap.compute_fft(sig, norm=True)

Pynapple returns a pandas DataFrame.

print(fft)
                          0
0.000   -0.002637+0.000000j
0.005   -0.003410-0.005622j
0.010    0.001544+0.009278j
0.015   -0.000156-0.000644j
0.020    0.000631+0.003293j
...                     ...
999.975  0.001622+0.007456j
999.980  0.006611+0.002686j
999.985 -0.001544-0.002222j
999.990 -0.001423+0.002182j
999.995  0.011593+0.011129j

[200000 rows x 1 columns]

It is then easy to plot it.

plt.figure()
plt.plot(np.abs(fft))
plt.xlabel("Frequency (Hz)")
plt.ylabel("Amplitude")
Text(0, 0.5, 'Amplitude')
../_images/5153b23fd88ab02d15b090278542613289efdd811240b505a189b76d4027b64f.png

Note that the output of the FFT is truncated to positive frequencies. To get positive and negative frequencies, you can set full_range=True. By default, the function returns the frequencies up to the Nyquist frequency. Let’s zoom on the first 20 Hz.

plt.figure()
plt.plot(np.abs(fft))
plt.xlabel("Frequency (Hz)")
plt.ylabel("Amplitude")
plt.xlim(0, 20)
(0.0, 20.0)
../_images/d221ee0368c7c4c398929af40b6c25844849988cf93478dafff85a5661e1f744.png

We find the two frequencies 2 and 10 Hz.

By default, pynapple assumes a constant sampling rate and a single epoch. For example, computing the FFT over more than 1 epoch will raise an error.

double_ep = nap.IntervalSet([0, 50], [20, 100])

try:
    nap.compute_fft(sig, ep=double_ep)
except ValueError as e:
    print(e)
Given epoch (or signal time_support) must have length 1

Computing power spectral density (PSD)#

Power spectral density can be returned through the function compute_power_spectral_density. Contrary to compute_fft, the output is real-valued.

psd = nap.compute_power_spectral_density(sig, fs=Fs)

The output is also a pandas DataFrame.

plt.figure()
plt.plot(psd)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Power/Frequency")
plt.xlim(0, 20)
(0.0, 20.0)
../_images/c8703de6cb4ae798ff3ae01bffb90fb7e5eba8f91417a3409f5525d281a530b8.png

Computing mean PSD#

It is possible to compute an average PSD over multiple epochs with the function nap.compute_mean_power_spectral_density.

In this case, the argument interval_size determines the duration of each epochs upon which the FFT is computed. If not epochs is passed, the function will split the time_support.

In this case, the FFT will be computed over epochs of 20 seconds.

mean_psd = nap.compute_mean_power_spectral_density(sig, interval_size=20.0, fs=Fs)

Let’s compare mean_psd to psd. In both cases, the output is normalized and converted to dB.

Hide code cell source
plt.figure()
plt.plot(10*np.log10(psd), label='PSD')
plt.plot(10*np.log10(mean_psd), label='Mean PSD (20s)')

plt.ylabel("Power/Frequency (dB/Hz)")
plt.legend()
plt.xlim(0, 20)
plt.xlabel("Frequency (Hz)")
Text(0.5, 0, 'Frequency (Hz)')
../_images/799e4c92432aac0dd3360b6f686df8b1eff582074a0190e87fb5f650aa1adb19.png

As we can see, nap.compute_mean_power_spectral_density was able to smooth out the noise.