Power spectral density

Show code cell source 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])+np.random.normal(0, 3, len(t)),
    time_support = nap.IntervalSet(0, 200)
    )

Show code cell source 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)')

Computing power spectral density (PSD)

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

psd = nap.compute_power_spectral_density(sig, norm=True)

Pynapple returns a pandas DataFrame.

print(psd)

                          0
0.000    0.000861+0.000000j
0.005   -0.003518+0.000328j
0.010    0.002092-0.000953j
0.015   -0.005350-0.002686j
0.020   -0.003940+0.000515j
...                     ...
999.975  0.000725-0.003469j
999.980  0.002271-0.003696j
999.985  0.000349+0.001222j
999.990 -0.002320+0.005219j
999.995 -0.000531-0.006738j

[200000 rows x 1 columns]

It is then easy to plot it.

plt.figure()
plt.plot(np.abs(psd))
plt.xlabel("Frequency (Hz)")
plt.ylabel("Amplitude")

Text(0, 0.5, 'Amplitude')

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(psd))
plt.xlabel("Frequency (Hz)")
plt.ylabel("Amplitude")
plt.xlim(0, 20)

(0.0, 20.0)

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_power_spectral_density(sig, ep=double_ep)
except ValueError as e:
    print(e)

Given epoch (or signal time_support) must have length 1

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 10 seconds.

mean_psd = nap.compute_mean_power_spectral_density(sig, interval_size=20.0, norm=True)

Let’s compare mean_psd to psd. In both cases, the ouput is normalized.

Show code cell source Hide code cell source

plt.figure()
plt.plot(np.abs(psd), label='PSD')
plt.plot(np.abs(mean_psd), label='Mean PSD (10s)')
plt.xlabel("Frequency (Hz)")
plt.ylabel("Amplitude")
plt.legend()
plt.xlim(0, 15)

(0.0, 15.0)

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