Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bumping v0.8.3 #403

Merged
merged 6 commits into from
Jan 25, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
20 changes: 10 additions & 10 deletions doc/examples/tutorial_wavelet_decomposition.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)
```
Expand All @@ -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()
```

Expand All @@ -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.

Expand All @@ -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
)
```

Expand All @@ -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()
```
Expand Down
4 changes: 4 additions & 0 deletions doc/releases.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
54 changes: 38 additions & 16 deletions doc/user_guide/10_power_spectral_density.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)
)
```
Expand All @@ -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")
```
Expand All @@ -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)
Expand All @@ -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
------------------

Expand All @@ -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.
22 changes: 11 additions & 11 deletions doc/user_guide/12_filtering.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)
```

Expand Down Expand Up @@ -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.

Expand Down
2 changes: 1 addition & 1 deletion pynapple/__init__.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
__version__ = "0.8.2"
__version__ = "0.8.3"
from .core import (
IntervalSet,
Ts,
Expand Down
6 changes: 5 additions & 1 deletion pynapple/process/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
Loading
Loading