Power Spectrum of CBGTC spike times

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

Bases: object

Computes the power spectra (Welch’s method using scipy.signal.welch) of all the neurons with given spike times

Methods

Argument

compute_for_spike()
  • spiketimes_set: see get_spiketimes_superset

  • sampling_rate [OPTIONAL]: 1000/dt = 10000 Hz [default]

  • window [OPTIONAL]: Tuple (0, 10) seconds [default]

  • neurons [OPTIONAL]: “all” [default] or a scalar or list: range(a, b) or [1, 4, 5, 9]

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

compute_for_rate()

  • mu_rate_array: see mean_rate()

  • binsz: integer or float

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

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

Use Cases

1. Pre-requisites

1.1. Import Modules

from analyseur.cbgtc.loader import LoadSpikeTimes
from analyseur.cbgtc.stats.psd import PowerSpectrum

1.2. Load file and get spike times

loadST = LoadSpikeTimes("spikes_GPi.csv")
spiketimes_set = loadST.get_spiketimes_superset()

2. Cases

2.1. Compute power spectral density (for all neurons)

B = PowerSpectrum.compute(spiketimes_set)

2.2. Compute power spectral density for chosen neurons with desired frequency resolution

B = PowerSpectrum.compute(spiketimes_set, neurons=range(30, 120), resolution=5)

Power spectral density for neurons 30 to 120 with the desired frequency resolution of 5 Hz.

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 given neurons.

\[P_r(f_k) = \frac{1}{N f_s} \left| \sum_{t=0}^{N-1} r(t)\, e^{-i 2\pi k t / N} \right|^2\]

is the squared magnitude of the Fourier transform of the population firing rate signal \(r(t)\) with number of samples \(N\) and sampling frequency is \(f_s = 1 / \Delta t\).

Parameters:

  • mu_rate_array – array of average firing rates for all/multiple neurons using mean_rate()

  • 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)

where the constant of proportionality is the number of points per segment nperseg.

Note that the population rate signal is

\[r(t) = \frac{1}{N}\sum_{i=1}^N s_i(t)\]

where \(i\)-th neuron has spike train \(s_i(t)\). Therefore, using the Fourier transform operator \(\mathcal{F}\) the power spectrum can be re-written as

\[P_r(f) = \left|\mathcal{F}\{r(t)\}\right|^2\]

Population synchrony at oscillatory frequency \(f_0\)

oscillation \(r(t) = A \cdot sin(2\pi f_0 t)\)

PSD peak, \(P_r(f_0)\)
classmethod compute_for_spike(spiketimes_set, sampling_rate=None, window=None, neurons=None, resolution=None)[source]

Returns the power spectral density (or power spectrum) of spiking times from all given neurons.

\[P_i(f) = \frac{1}{K} \sum_{k=1}^K \left[\frac{1}{L}\left|\sum_{n=0}^{L-1}s_i(n)\cdot w(n) \cdot e^{-i2\pi f n}\right|^2\right]\]

is the Welch’s estimate of the power spectral density \(P_i(f)\) with \(w(n)\) window, \(K\) number of segments, \(L\) nperseg, and the binary spike train \(s_i(t_k)\) (see get_binary_spiketrains()).

Parameters:

spiketimes_set – Dictionary returned using get_spiketimes_superset()

or using get_spiketimes_subset()

Parameters:

  • sampling_rate – 1000/dt = 10000 Hz [default]; sampling_rate ∊ (0, 10000)

  • window – Tuple in the form (start_time, end_time); (0, 10) [default]

  • neurons –

    “all” [default] or scalar or range(a, b) or list of neuron ids like [2, 3, 6, 7]

    • ”all” means subset = superset

    • N (a scalar) means subset of first N neurons in the superset

    • range(a, b) or [2, 3, 6, 7] means subset of selected neurons

  • 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)

  • list of spike trains

  • list of neuron id’s

  • array of time

Notes on resolution

resolution or desired frequency resolution is the smallest difference between two frequencies that can be distinguished in the power spectrum.

The sampling rate is proportional to the desired frequency resolution.
\[ \begin{align}\begin{aligned}\text{sampling\_rate} &\propto \text{resolution} \\\text{sampling\_rate} &= \text{nperseg} \times \text{resolution}\end{aligned}\end{align} \]

where the constant of proportionality is the number of points per segment nperseg.

Analysis Pitfalls

\(P_i(f)\) returns one PSD per neuron

plotting PSD for just one neuron will mostly be Poisson noise

the meaningful quantity is \(P_\text{pop}(f) = \frac{1}{N}\sum_{i=1}^{N}P_i(f)\)

\(P_\text{pop}(f)\) is used in plot_aggstat()