Spectral analysis & continous data processing (spynal.spectra)

Spectral analysis, signal processing, and continuous (LFP/EEG) data preprocessing

This module is the base module for the spectra module. It contains all the high-level API functions. Lower-level functions can be found in the other spectra modules.

Overview

Functionality for computing frequency spectra as well as time-frequency (spectrogram) transforms.

Options to compute spectral analysis using multitaper, wavelet, band-pass filtering, or spectral burst analysis methods.

Options to return full complex spectral data, spectral power, phase, real or imaginary part, etc.

Also includes functions for preprocessing, postprocessing, plotting of continuous/spectral data.

Most functions perform operations in a mass-univariate manner. This means that rather than embedding function calls in for loops over channels, trials, etc., like this:

for channel in channels:
    for trial in trials:
        results[trial,channel] = compute_something(data[trial,channel])

You can instead execute a single call on ALL the data, labeling the relevant axis for the computation (usually time here), and it will run in parallel (vectorized) across all channels, trials, etc. in the data, like this:

results = compute_something(data, axis)

Function list

General spectral analysis

  • spectrum : Frequency spectrum of data

  • spectrogram : Time-frequency spectrogram of data

  • power_spectrum : Power spectrum of data

  • power_spectrogram : Power of time-frequency transform

  • phase_spectrogram : Phase of time-frequency transform

Wavelet spectral analysis (spectra.wavelet)

  • wavelet_spectrum : Wavelet-based frequency spectrum

  • wavelet_spectrogram : Time-frequency continuous wavelet transform

  • compute_wavelets : Compute wavelets for use in wavelet spectral analysis

  • wavelet_bandwidth : Compute time,frequency bandwidths for set of wavelets

  • wavelet_edge_extent : Compute extent of edge effects for set of wavelets

Multitaper spectral analysis (spectra.multitaper)

  • multitaper_spectrum : Multitaper (DPSS) frequency spectrum

  • multitaper_spectrogram : Multitaper (DPSS) time-frequency spectrogram

  • compute_tapers : Compute DPSS tapers for use in multitaper spectral analysis

Bandpass-filtering spectral analysis (spectra.bandfilter)

  • bandfilter_spectrum : Band-filtered frequency spectrum

  • bandfilter_spectrogram : Band-filtered, Hilbert-transformed time-frequency of data

  • set_filter_params : Set filter coefficients for use in band-filtered analysis

Other spectral analyses

  • itpc : Intertrial phase clustering (analysis of phase locking to trial events)

  • burst_analysis : Compute oscillatory burst analysis of Lundqvist et al 2016

Preprocessing (spectra.preprocess)

  • cut_trials : Cut LFPs/continuous data into trial segments

  • realign_data : Realign LFPs/continuous data to new within-trial event

  • remove_dc : Remove constant DC component of signals

  • remove_evoked : Remove phase-locked evoked potentials from signals

Postprocesssing (spectra.postprocess)

  • pool_freq_bands : Average spectral data within set of frequency bands

  • pool_time_epochs : Average spectral data within set of time epochs

  • one_over_f_norm : Normalize to correct for 1/f distribution of spectral power

Utilities (spectra.utils)

  • get_freq_sampling : Frequency sampling vector for a given FFT-based computation

  • complex_to_spec_type : Convert complex Fourier transform output to power/phase/real/imag/etc.

  • one_sided_to_two_sided : Convert 1-sided Fourier transform output to 2-sided equivalent

  • simulate_oscillation : Generates simulated oscillation-in-noise data

Plotting

  • plot_spectrum : Plot frequency spectrum as a line plot, handling freq axis properly

  • plot_spectrogram : Plot time-frequency spectrogram as a heatmap plot

Dependencies

  • pyfftw : Python wrapper around FFTW, the speedy FFT library

Function reference

spectrum(data, smp_rate, axis=0, method='multitaper', data_type='lfp', spec_type='complex', removeDC=True, **kwargs)

Compute frequency spectrum of data using given method

Parameters:

  • data (ndarray,shape=(...,n_samples,...)) – Data to compute spectral analysis of. Arbitrary shape; spectral analysis is computed along axis.

  • smp_rate (scalar) – Data sampling rate (Hz)

  • axis (int, default: 0 (1st axis)) – Axis of data to do spectral analysis on (usually time dimension).

  • method ({'multitaper','wavelet','bandfilter'}, default: 'multitaper') –

    Specific spectral analysis method to use:

  • data_type ({'lfp','spike'}, default: 'lfp') – Type of signal in data

  • spec_type ({'complex','power','phase','real','imag'}, default: 'complex') – Type of spectral signal to return. See complex_to_spec_type() for details.

  • removeDC (bool, default: True) – If True, subtracts off mean DC component across axis, making signals zero-mean before spectral analysis.

  • **kwargs – All other kwargs passed directly to method-specific spectrum function

Returns:

  • spec (ndarray, shape=(…,n_freqs,…), dtype=complex or float.) – Frequency spectrum of given type computed with given method. Frequency axis is always inserted in place of axis. Note: ‘multitaper’ method will return with additional taper axis inserted after just after axis if keep_tapers is True. dtype is complex if spec_type is ‘complex’, float otherwise.

  • freqs (ndarray, shape=(n_freqs,) or (n_freqbands,2)) – For method == ‘bandfilter’: List of (low,high) cut frequencies (Hz) used to generate spec, shape=(n_freqbands,2) For other methods: List of frequencies in spec (Hz), shape=(n_freqs,)

spectrogram(data, smp_rate, axis=0, method='wavelet', data_type='lfp', spec_type='complex', removeDC=True, **kwargs)

Compute time-frequency transform of data using given method

Parameters:

  • data (ndarray, shape=(...,n_samples,...)) – Data to compute spectral analysis of. Arbitrary shape; spectral analysis is computed along axis.

  • smp_rate (scalar) – Data sampling rate (Hz)

  • axis (int, default: 0 (1st axis)) – Axis of data to do spectral analysis on (usually time dimension).

  • method ({'multitaper','wavelet','bandfilter','burst'}, default: 'wavelet') –

    Specific spectral analysis method to use:

  • data_type ({'lfp','spike'}, default: 'lfp') – Type of signal in data

  • spec_type ({'complex','power','phase','real','imag'}, default: 'complex') – Type of spectral signal to return. See complex_to_spec_type() for details.

  • removeDC (bool, default: True) – If True, subtracts off mean DC component across axis, making signals zero-mean before spectral analysis.

  • **kwargs – All other kwargs passed directly to method-specific spectrogram function

Returns:

  • spec (ndarray, shape=(…,n_freqs,n_timepts,…), dtype=complex or float.) – Time-frequency spectrogram of given type computed with given method. Frequency axis is always inserted just before time axis. Note: ‘multitaper’ method will return with additional taper axis inserted between freq and time axes if keep_tapers is True. dtype is complex if spec_type is ‘complex’, float otherwise.

  • freqs (ndarray, shape=(n_freqs,) or (n_freqbands,2)) – For method == ‘bandfilter’: List of (low,high) cut frequencies (Hz) used to generate spec, shape=(n_freqbands,2) For other methods: List of frequencies in spec (Hz), shape=(n_freqs,)

  • timepts (ndarray, shape=(n_timepts,)) – List of time points / time window centers in spec (in s, referenced to start of data)

power_spectrum(data, smp_rate, axis=0, method='multitaper', **kwargs)

Convenience wrapper around spectrum() to compute power spectrum of data with given method

See spectrum() for details

power_spectrogram(data, smp_rate, axis=0, method='wavelet', **kwargs)

Convenience wrapper around spectrogram() to compute time-frequency power with given method

See spectrogram() for details

phase_spectrogram(data, smp_rate, axis=0, method='wavelet', **kwargs)

Convenience wrapper around spectrogram() to compute phase of time-frequency transform

See spectrogram() for details

itpc(data, smp_rate, axis=0, method='wavelet', itpc_method='PLV', trial_axis=None, **kwargs)

Intertrial phase clustering (ITPC) measures frequency-specific phase locking of continuous neural activity (eg LFPs) to trial events. A spectral analog (roughly) of evoked potentials.

Complex time-frequency representation is first computed (using some method), then the complex vector mean is computed across trials, and it’s magnitude is returned as ITPC.

aka “intertrial coherence”, “intertrial phase-locking value/factor”

Parameters:

  • data (ndarray, shape=(...,n_timepts,...,n_trials,...)) – Data to compute ITPC of, aligned to some within-trial event. Arbitrary shape; spectral analysis is computed along axis (usually time), and ITPC is computed along trial_axis.

  • smp_rate (scalar) – Data sampling rate (Hz)

  • axis (int, default: 0 (1st axis)) – Axis of data to do spectral analysis on (usually time dimension).

  • method ({'multitaper','wavelet','bandfilter'}, default: 'wavelet') –

    Specific underlying spectral analysis method to use:

  • itpc_method ({'PLV','Z','PPC'}, default: 'PLV') –

    Method to use for computing intertrial phase clustering:

    • ’PLV’Phase locking value (length of cross-trial complex vector mean)

      Standard/traditional measure of ITPC, but is positively biased for small n, and is biased > 0 even for no phase clustering.

    • ’Z’ : Rayleigh’s Z normalization of PLV (Z = n*PLV**2). Reduces small-n bias.

    • ’PPC’Pairwise Phase Consistency normalization (PPC = (n*PLV**2 - 1)/(n - 1))

      of PLV. Debiased and has expected value 0 for no clustering.

  • **kwargs – All other kwargs passed directly to method-specific spectrogram function

Returns:

  • ITPC (ndarray, shape=(…,n_freqs,n_timepts,…)) – Time-frequency spectrogram representation of intertrial phase clustering. Frequency axis is always inserted just before time axis, and trial axis is removed, but otherwise shape is same as input data.

  • freqs (ndarray, shape=(n_freqs,) or (n_freqbands,2)) – For method == ‘bandfilter’: List of (low,high) cut frequencies (Hz) used to generate ITPC, shape=(n_freqbands,2) For other methods: List of frequencies in spec (Hz), shape=(n_freqs,)

  • timepts (ndarray, shape=(n_timepts,)) – List of time points / time window centers in ITPC (in s, referenced to start of data)

References

Cohen “Analyzing Neural Time Series Data” http://dx.doi.org/10.7551/mitpress/9609.001.0001 Ch. 19

intertrial_phase_clustering(data, smp_rate, axis=0, method='wavelet', itpc_method='PLV', trial_axis=None, **kwargs)

Alias of itpc(). See there for details.

burst_analysis(data, smp_rate, axis=0, trial_axis=-1, threshold=2, min_cycles=3, method='wavelet', spec_type='power', freq_exp=None, bands=((20, 35), (40, 65), (55, 90), (70, 100)), window=None, timepts=None, **kwargs)

Oscillatory burst analysis of Lundqvist et al 2016.

Computes oscillatory power, z-scores within each frequency band, thresholds at given z threshold, labels as burst “ON” times timepoints > threshold for at least min_cycles duration.

To compute burst rate, simply take mean of computed bursts across trial axis.

Default argument values approximate analysis as implemented in Lundqvist 2016.

Parameters:

  • data (ndarray, shape=(...,n_samples,...)) – Data to compute spectral analysis of. Arbitrary shape; spectral analysis is computed along axis.

  • smp_rate (scalar) – Data sampling rate (Hz)

  • smp_rate – Data sampling rate (Hz)

  • axis (int, default: 0 (1st axis)) – Axis of data to do spectral analysis on (usually time dimension).

  • trial_axis (int, default: -1 (last axis of data)) – Axis of data corresponding to trials/observations.

  • threshold (scalar, default: 2 (2 SDs above mean)) – Threshold power level for detecting bursts, given in SDs above the mean

  • min_cycles (scalar, default: 3) – Minimal length of contiguous above-threshold period to be counted as a burst, given in number of oscillatory cycles at each frequency (or band center).

  • method ({'multitaper','wavelet','bandfilter'}, default: 'wavelet') –

    Specific spectral analysis method to use:

    Note: In the original paper, multitaper was used, but all three methods were claimed to produced similar results.

  • spec_type ({'power','magnitude'}, default: 'power') –

    Type of spectral signal to compute:

    • ’power’ : Spectral power, ie square of signal envelope

    • ’magnitude’ : Square root of power, ie signal envelope

  • freq_exp (float, default: None (do no frequency normalization)) – This can be used to normalize out 1/f^a effects in power before band-pooling and burst detection). This gives the exponent on the frequency (‘a’ in 1/f^a). Set = 1 to norm by 1/f. Set = None for no normalization.

  • bands (array-like, shape=(n_freqbands,2), default: ((20,35),(40,65),(55,90),(70,100))) – List of (low,high) cut frequencies for each band to compute bursts within. Set low cut = 0 for low-pass, set high cut >= smp_rate/2 for high-pass, otherwise assumes band-pass. Default samples ~ beta, low/med/high gamma. Set = None to compute bursts at each frequency in spectral transform.

  • window (array-like, shape=(2,), default: None (compute over entire data time range)) – Optional (start,end) of time window to compute mean,SD for burst amplitude threshold within (in same units as timepts).

  • timepts (array-like, shape=(n_timepts,), default: None) – Time sampling vector for data (usually in s). Necessary if window is set, but unused otherwise.

  • **kwargs – Any other kwargs passed directly to spectrogram() function

Returns:

  • bursts (ndarray, shape=(…,n_freq[band]s,n_timepts_out,…), dtype=bool) – Binary array labelling timepoints within detected bursts in each trial and frequency band. Same shape as data, but with frequency axis prepended immediately before time axis.

  • freqs (ndarray, shape=(n_freq[band],)) – List of center frequencies of bands in bursts

  • timepts (ndarray=(n_timepts_out,)) – List of time points / time window centers in bursts (in s, referenced to start of data).

References

Lundqvist et al (2016) Neuron https://doi.org/10.1016/j.neuron.2016.02.028 Lundqvist et al (2018) Nature Comms https://doi.org/10.1038/s41467-017-02791-8

plot_spectrum(freqs, data, ax=None, ylim=None, color=None, **kwargs)

Plot frequency spectrum as a line plot.

Parameters:

  • freqs (array-like, shape=(n_freqs,)) – Frequency sampling (x-axis) vector for data (Hz). May be linearly or logarithmically sampled; we handle appropriately.

  • data (ndarray, shape=(n_freqs,)) – Frequency spectrum data to plot (y-axis)

  • ax (Pyplot Axis object, default: plt.gca()) – Axis to plot into

  • ylim (array-like, shape=(2,), Default: (data.min(),data.max()) +/- 5%) – Plot y-axis limits: (min,max)

  • color (Color spec, default: <Matplotlib default plot color>) – Color to plot in

  • **kwargs – Any additional keyword args are interpreted as parameters of plt.axes() (settable Axes object attributes) or plt.plot() (Line2D object attributes), and passsed to the proper function.

Returns:

  • lines (List of Line2D objects) – Output of plt.plot()

  • ax (Axis object) – Axis plotted into

plot_spectrogram(timepts, freqs, data, ax=None, clim=None, cmap='viridis', **kwargs)

Plot time-frequency spectrogram as a heatmap plot.

Parameters:

  • timepts (array-like, shape=(n_timepts,)) – Time sampling (x-axis) vector for data

  • freqs (array-like, shape=(n_freqs,)) – Frequency sampling (y-axis) vector for data (Hz). May be linearly or logarithmically sampled; we handle appropriately.

  • data (ndarray, shape=(n_freqs,n_timepts)) – Time-frequency (spectrogam) data to plot on color axis

  • ax (Pyplot Axis object, default: plt.gca()) – Axis to plot into

  • clim (array-like, shape=(2,), default: (data.min(),data.max())) – Color axis limits: (min,max)

  • cmap (str | Colormap object. default: 'viridis' (linear dark-blue to yellow colormap)) – Colormap to plot heatmap in, given either as name of matplotlib colormap or custom matplotlib.colors.Colormap object instance.

  • **kwargs – Any additional keyword args are interpreted as parameters of plt.axes() (settable Axes object attributes) or plt.imshow() (AxesImage object attributes).

Returns:

  • img (AxesImage object) – Output of ax.imshow(). Allows access to image properties.

  • ax (Axis object) – Axis plotted into.