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/77e810daf274ec230a283b02584160620fb0c02c3c073fbf29a85e5ef5594d31.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.002519+0.000000j
0.005    0.008462-0.005830j
0.010   -0.005360-0.001660j
0.015    0.001949+0.003159j
0.020    0.010565-0.009647j
...                     ...
999.975 -0.005571+0.002369j
999.980  0.006477-0.010829j
999.985 -0.001354+0.001348j
999.990 -0.004499+0.003365j
999.995  0.000163+0.003419j

[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/e8d4bc3c4521a03f1ebea58f1719a286433be3d5e627c557c4e746d428e9e6f3.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/9fe967de94e304d2b31d23a93d1e21d0fac1c76a5e7a53a36d3540ad6ac16e73.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/b7ef9eb3ef750a47b7271a0004417895b93f6ba976afc2cc715e66ab5c1e1508.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/0e212c97e6c105902f03198d2a31bb72e222cbb499f7b49db7a9f657882550f0.png

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