Merge pull request #403 from pynapple-org/dev

Bumping v0.8.3

doc/examples/tutorial_wavelet_decomposition.md

@@ -151,7 +151,7 @@ Let's take the Fourier transform of our data to get an initial insight into the
 jupyter:
   outputs_hidden: false
 ---
-power = nap.compute_power_spectral_density(eeg, fs=FS, ep=wake_ep, norm=True)
+power = nap.compute_power_spectral_density(eeg, fs=FS, ep=wake_ep)
 
 print(power)
 ```
@@ -169,14 +169,14 @@ jupyter:
 ---
 fig, ax = plt.subplots(1, constrained_layout=True, figsize=(10, 4))
 ax.plot(
-    np.abs(power[(power.index >= 1.0) & (power.index <= 100)]),
+    power[(power.index >= 1.0) & (power.index <= 100)],
     alpha=0.5,
     label="LFP Frequency Power",
 )
 ax.axvspan(6, 10, color="red", alpha=0.1)
 ax.set_xlabel("Freq (Hz)")
-ax.set_ylabel("Frequency Power")
-ax.set_title("LFP Fourier Decomposition")
+ax.set_ylabel("Power/Frequency ")
+ax.set_title("LFP Power spectral density")
 ax.legend()
 ```
 
@@ -200,7 +200,7 @@ rest_ep = wake_ep.set_diff(run_ep)
 
 The function `nap.compute_power_spectral_density` takes signal with a single epoch to avoid artefacts between epochs jumps.
 
-To compare `run_ep` with `rest_ep`, we can use the function `nap.compute_mean_power_spectral_density` which avearge the FFT over multiple epochs of same duration. The parameter `interval_size` controls the duration of those epochs.
+To compare `run_ep` with `rest_ep`, we can use the function `nap.compute_mean_power_spectral_density` which average the FFT over multiple epochs of same duration. The parameter `interval_size` controls the duration of those epochs.
 
 In this case, `interval_size` is equal to 1.5 seconds.
 
@@ -211,10 +211,10 @@ jupyter:
   outputs_hidden: false
 ---
 power_run = nap.compute_mean_power_spectral_density(
-    eeg, 1.5, fs=FS, ep=run_ep, norm=True
+    eeg, 1.5, fs=FS, ep=run_ep
 )
 power_rest = nap.compute_mean_power_spectral_density(
-    eeg, 1.5, fs=FS, ep=rest_ep, norm=True
+    eeg, 1.5, fs=FS, ep=rest_ep
 )
 ```
 
@@ -228,20 +228,20 @@ jupyter:
 ---
 fig, ax = plt.subplots(1, constrained_layout=True, figsize=(10, 4))
 ax.plot(
-    np.abs(power_run[(power_run.index >= 3.0) & (power_run.index <= 30)]),
+    power_run[(power_run.index >= 3.0) & (power_run.index <= 30)],
     alpha=1,
     label="Run",
     linewidth=2,
 )
 ax.plot(
-    np.abs(power_rest[(power_rest.index >= 3.0) & (power_rest.index <= 30)]),
+    power_rest[(power_rest.index >= 3.0) & (power_rest.index <= 30)],
     alpha=1,
     label="Rest",
     linewidth=2,
 )
 ax.axvspan(6, 10, color="red", alpha=0.1)
 ax.set_xlabel("Freq (Hz)")
-ax.set_ylabel("Frequency Power")
+ax.set_ylabel("Power/Frequency")
 ax.set_title("LFP Fourier Decomposition")
 ax.legend()
 ```

doc/releases.md

@@ -18,6 +18,10 @@ of the Flatiron institute.
 
 ## Releases
 
+### 0.8.3 (2025-01-24)
+
+- `compute_mean_power_spectral_density` computes the mean periodogram.
+
 ### 0.8.2 (2025-01-22)
 
 - `compute_power_spectral_density` now computes the periodogram, where previously it was only computing the FFT

doc/user_guide/10_power_spectral_density.md

@@ -40,7 +40,7 @@ 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)),
+    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)
     )
 ```
@@ -54,29 +54,29 @@ plt.title("Signal")
 plt.xlabel("Time (s)")
 ```
 
-Computing power spectral density (PSD)
+Computing FFT
 --------------------------------------
 
-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.
+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.
 
 
 ```{code-cell} ipython3
-psd = nap.compute_power_spectral_density(sig, norm=True)
+fft = nap.compute_fft(sig, norm=True)
 ```
 
 Pynapple returns a pandas DataFrame.
 
 
 ```{code-cell} ipython3
-print(psd)
+print(fft)
 ```
 
 It is then easy to plot it.
 
 
 ```{code-cell} ipython3
 plt.figure()
-plt.plot(np.abs(psd))
+plt.plot(np.abs(fft))
 plt.xlabel("Frequency (Hz)")
 plt.ylabel("Amplitude")
 ```
@@ -88,7 +88,7 @@ Let's zoom on the first 20 Hz.
 
 ```{code-cell} ipython3
 plt.figure()
-plt.plot(np.abs(psd))
+plt.plot(np.abs(fft))
 plt.xlabel("Frequency (Hz)")
 plt.ylabel("Amplitude")
 plt.xlim(0, 20)
@@ -103,11 +103,31 @@ By default, pynapple assumes a constant sampling rate and a single epoch. For ex
 double_ep = nap.IntervalSet([0, 50], [20, 100])
 
 try:
-    nap.compute_power_spectral_density(sig, ep=double_ep)
+    nap.compute_fft(sig, ep=double_ep)
 except ValueError as e:
     print(e)
 ```
 
+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.
+
+```{code-cell} ipython3
+psd = nap.compute_power_spectral_density(sig, fs=Fs)
+```
+
+The output is also a pandas DataFrame.
+
+```{code-cell} ipython3
+plt.figure()
+plt.plot(psd)
+plt.xlabel("Frequency (Hz)")
+plt.ylabel("Power/Frequency")
+plt.xlim(0, 20)
+```
+
 Computing mean PSD
 ------------------
 
@@ -116,26 +136,28 @@ It is possible to compute an average PSD over multiple epochs with the function
 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.
+In this case, the FFT will be computed over epochs of 20 seconds.
 
 
 ```{code-cell} ipython3
-mean_psd = nap.compute_mean_power_spectral_density(sig, interval_size=20.0, norm=True)
+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 ouput is normalized.
+Let's compare `mean_psd` to `psd`. In both cases, the output is normalized and converted to dB.
 
 
 ```{code-cell} ipython3
 :tags: [hide-input]
 
 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.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, 15)
+plt.xlim(0, 20)
+plt.xlabel("Frequency (Hz)")
+
 ```
 
 As we can see, `nap.compute_mean_power_spectral_density` was able to smooth out the noise.

doc/user_guide/12_filtering.md

@@ -72,15 +72,15 @@ We can compute the Fourier transform of `sig` to verify that all the frequencies
 
 
 ```{code-cell} ipython3
-psd = nap.compute_power_spectral_density(sig, fs, norm=True)
+psd = nap.compute_power_spectral_density(sig, fs)
 ```
 ```{code-cell} ipython3
 :tags: [hide-input]
 
 fig = plt.figure(figsize = (15, 5))
-plt.plot(np.abs(psd))
+plt.plot(psd)
 plt.xlabel("Frequency (Hz)")
-plt.ylabel("Amplitude")
+plt.ylabel("Power/Frequency")
 plt.xlim(0, 100)
 ```
 
@@ -146,29 +146,29 @@ Let's see what frequencies remain;
 
 
 ```{code-cell} ipython3
-psd_butter = nap.compute_power_spectral_density(sig_butter, fs, norm=True)
-psd_sinc = nap.compute_power_spectral_density(sig_sinc, fs, norm=True)
+psd_butter = nap.compute_power_spectral_density(sig_butter, fs)
+psd_sinc = nap.compute_power_spectral_density(sig_sinc, fs)
 ```
 
 ```{code-cell} ipython3
 :tags: [hide-input]
 
 fig = plt.figure(figsize = (10, 5))
-plt.plot(np.abs(psd_butter), label = "Butterworth filter")
-plt.plot(np.abs(psd_sinc), label = "Windowed-sinc convolution")
+plt.plot(psd_butter, label = "Butterworth filter")
+plt.plot(psd_sinc, label = "Windowed-sinc convolution")
 plt.legend()
 plt.xlabel("Frequency (Hz)")
-plt.ylabel("Amplitude")
+plt.ylabel("Power/Frequency")
 plt.xlim(0, 70)
 ```
 
 ***
-Inspecting frequency fesponses of a filter
+Inspecting frequency responses of a filter
 ------------------------------------------
 
-We can inspect the frequency response of a filter by plotting its power spectral density (PSD).
+We can inspect the frequency response of a filter by plotting its FFT.
 To do this, we can use the `get_filter_frequency_response` function, which returns a pandas Series with the frequencies
-as the index and the PSD as values.
+as the index and the amplitude as values.
 
 Let's extract the frequency response of a Butterworth filter and a sinc low-pass filter.
 

pynapple/__init__.py

@@ -1,4 +1,4 @@
-__version__ = "0.8.2"
+__version__ = "0.8.3"
 from .core import (
     IntervalSet,
     Ts,

pynapple/process/__init__.py

@@ -22,7 +22,11 @@
     shift_timestamps,
     shuffle_ts_intervals,
 )
-from .spectrum import compute_fft, compute_mean_fft, compute_power_spectral_density
+from .spectrum import (
+    compute_fft,
+    compute_mean_power_spectral_density,
+    compute_power_spectral_density,
+)
 from .tuning_curves import (
     compute_1d_mutual_info,
     compute_1d_tuning_curves,