Power Spectrum of rBCBG firing rates

class analyseur.rbcbg.stats.psd.PowerSpectrum[source]

Bases: object

Computes the power spectra (Welch’s method using scipy.signal.welch) of the firing rates for all the neurons.

Methods

Argument

compute_spectrogram()

  • mu_rate_array: see get_rates()

  • resolution [OPTIONAL]: ~ 9.76 Hz = sampling_rate/1024 [default]

compute_for_rate()

  • mu_rate_array: see get_rates()

  • binsz: integer or float

  • method: “welch” or “fft” or “fft-mag”

  • resolution [OPTIONAL]: ~ 9.76 Hz = sampling_rate/1024 [default]

  • compute_spectrogram() computes a time-varying power spectrum, that is, power for each time window and hence the power spectrum evolving over time.

  • compute_for_rate() computes a single (global) power spectral density, that is, power averaged over the whole signal (no time dimension) and hence one summary spectrum for the entire signal.

Therefore, compute_spectrogram() retains the time structure while compute_for_rate() collapses time by averaging.

Use Cases

1. Pre-requisites

1.1. Import Modules

from analyseur.rbcbg.loader import LoadRates
from analyseur.rbcbg.stats.psd import PowerSpectrum

1.2. Load file and get the firing rates

loadFR = LoadRates("GPiSNr_model_9_percent_0.csv")
t_sec, rates_Hz = loadFR.get_rates()

2. Cases

2.1. Time-varying power spectrum for a desired resolution

freq_arr, time_arr, power_arr = PowerSpectrum.compute_spectrogram(rates_Hz, resolution=10)

done for resolution 10 Hz.

2.2. Global power spectral density

pow_spec = PowerSpectrum.compute_for_rate(mu_rate_array, binsz)

doe for default resolution and Welch method (default).

classmethod compute_for_rate(mu_rate_array, binsz, resolution=None, method=None)[source]

Returns the power spectral density (or power spectrum) of firing rate from all.

\[\hat{P}_{xx}(f) = \frac{1}{K} \sum_{m=1}^K S(m,f)\]

is the average of the spectrogram \(S(m,f)\) (see compute_spectrogram()) over time windows (Welch PSD). Thus, \(P(f)\) is averaged over \(m\) while \(S(m,f)\) depends on time.

Parameters:

  • mu_rate_array – array returned using get_rates()

  • binsz – integer or float

  • method – “welch” or “fft” or “fft-mag”

  • resolution – ~ 9.76 Hz = sampling_rate/1024 [default]

Returns:

a tuple in the following order

  • array of sample frequencies

  • power spectral density (or power spectrum)

From compute_spectrogram() substituting for \(S(m,f)\) the Welch PSD is expanded as

\[\hat{P}_{xx}(f) = \frac{1}{K} \sum_{m=1}^K \frac{1}{f_s||w||^2} \cdot \left|\sum_{n=0}^{N-1}x[n+m]\cdot w[n] \cdot e^{-j2\pi fn/f_s} \right|^2\]

where \(N\) is the segment length (nperseg), \(f_s\) is the sampling frequency (sample_rate) and \(||w||^2 = \sum_n w[n]^2\) is the normalized window.

Removing the windowing and segmentation we get the FFT PSD

\[P_{xx}(f) = \frac{1}{f_s N} \cdot \left|\sum_{n=0}^{N-1}x[n] \cdot e^{-j2\pi fn/f_s} \right|^2\]

This is the power spectrum from the single whole global window.

NOTE: This average (or single/global) power spectrum of the entire signal tells us the overall frequency content.

classmethod compute_spectrogram(rates_Hz, resolution=None)[source]

Compute spectrogram using STFT (see compute_stft())

Parameters:

  • rates_Hz – array returned using get_rates()

  • resolution – ~ 9.76 Hz = sampling_rate/1024 [default]

Returns:

a tuple in the following order

  • f: frequency vector

  • t: time vector for spectrogram

  • Sxx: spectrogram power (dB)

NOTE: This time-varying power spectrum tells us how frequencies evolve.

static compute_stft(rate_array, sample_rate, nperseg=256, noverlap=128)[source]

Computes short-time Fourier transform (STFT)

Parameters:

  • rate_array – array returned using get_rates()

  • sample_rate – float

  • nperseg – int, Length of each segment

  • noverlap – int, Overlaps between segments

Returns:

a tuple in the following order

  • f: frequency vector

  • t: time vector for spectrogram

  • Sxx: spectrogram power (dB)

Mathematically, the array of firing rates \(x[n]\) is segmented into overlapping windows and according to the window position \(m\) (time index) and angular frequency \(\omega\) it is Fourier transformed as

\[X(m,\omega) = \sum_{n=-\infty}^\infty x[n]\cdot w[n-m] \cdot e^{-j\omega n}\]

where the window function \(w[n]\) is given by the Hann window

\[w[n] = 0.5\left(1 - cos\left(\frac{2\pi n}{N-1}\right)\right)\]

where \(N\) is the segment length (nperseg).

Then the spectrogram is the squared magnitude of the STFT

\[ \begin{align}\begin{aligned}S(m,\omega) &= |X(m,\omega)|^2 \\S(m, f) &= \frac{1}{f_s||w||^2} \cdot \left|\sum_{n=0}^{N-1}x[n+m]\cdot w[n] \cdot e^{-j2\pi fn/f_s} \right|^2\end{aligned}\end{align} \]

where \(f_s\) is the sampling frequency (sample_rate) and \(||w||^2 = \sum_n w[n]^2\) is the normalized window.

Finally, for \(k=0,1,\ldots,N/2\) the frequency bins are

\[f_k = \frac{k}{N} \cdot f_s\]

the time bins are

\[t_m = \frac{m(N-\text{noverlap})}{f_s}\]

and the power is converted to decibels

\[S_\text{dB}(m,f) = 10 log_{10}(S(m,f) + \epsilon)\]

where \(\epsilon = 10^{-10}\) is a small constant to avoid \(log(0)\)

Note that the no overlap parameter noverlap when subtracted from the segment length \(N\) yields the hop size (step size) \(= N - \text{noverlap}\). It acts like a sliding window such that

noverlap = 0

[--------][--------][--------]

but with the presence of overlap

noverlap = 128

[--------]
    [--------]
        [--------]

resulting in more redundancy but better continuity and hence smoother time resolution.