General spectral analysis utilities (spynal.spectra.utils)

Utility functions for LFP/EEG/continuous data and spectral analysis

next_power_of_2(n)

Find next power of 2 (smallest power of 2 greater than n)

get_freq_sampling(smp_rate, n_fft, freq_range=None, pad=False, two_sided=False)

Return frequency sampling vector (axis) for a given FFT-based computation

Parameters:

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

• n_fft (scalar) – Number of samples (timepoints) in FFT output (excluding any padding due to pad)

• freq_range (array-like, shape=(2,) or scalar, default: all frequencies from FFT) – Range of frequencies to retain in output, either given as an explicit [low,high] range or just a scalar giving the highest frequency to return.

• pad (bool, default: False) – If True, n_fft is increased to length = the next power of 2

• two_sided (bool, default: False) – If True, return freqs for two-sided spectrum, including both positive and negative frequencies (which have same amplitude for all real signals). If False, only return positive frequencies, in range (0,smp_rate/2).

Returns:

• freqs (ndarray, shape=(n_freqs,)) – Frequency sampling vector (in Hz)

• freq_bool (ndarray, shape=(n_fft,), dtype=bool) – Boolean vector flagging frequencies in full FFT output to retain, given desired freq_range

get_freq_length(smp_rate, n_fft, freq_range=None, pad=False, two_sided=False)

Return length of frequency sampling vector (axis) for a given FFT-based computation.

Wrapper around get_freq_sampling() that simply returns the length of the list of sampled freqs. Input arguments are identical; see there for details.

Returns:

n_freqs – Length of frequency sampling vector for given FFT, which could be used to set axis length for array to hold FFT results.

Return type:

int

fft(data, n_fft=None, axis=0, fft_method=None)

Compute 1d discrete Fast Fourier transform along given array axis, using given method.

Thin wrapper around multiple low-level FFT implementations, which may be chosen based on performance and availability (in our hands, torch is fastest for virtually all tests).

Parameters:

• data (ndarray, shape=(...,n_samples,...)) – Data to compute 1d FFT on. Shape is arbitrary; FFT is performed along axis (usually corresponding to time) independently across any/all other array axes (which might be data channels, trials, etc.).

• n_fft (int, default: data.shape[axis]) – Length of the transformed axis of the output. If smaller than the length of the input, the input is cropped. If larger, the input is zero-padded with zeros. If not given, the length of the input along the axis specified by axis is used.

• axis (int, default: 0 (1st array axis)) – Axis over which to compute the FFT. If not given, the first axis is used.

• fft_method (str, default: 'torch' (if available)) –

  Which underlying low-level FFT implementation to use. Options:

  • ’torch’Torch’s FFT on the GPU, using :func:torch.fft.fft.

    Across a range of tests, this is by far the fastest method.

  • ’fftw’FFTW library’s optimized FFT, using :func:pyfftw.interfaces.scipy_fftpack.fft

    FFTW is faster than scipy/numpy for data length != power of 2 (but slower than torch).

  • ’numpy’ : Numpy’s FFT implementation, using :func:np.fft.fft

  • ’scipy’ : Scipy’s FFT implementation, using :func:scipy.fft.fft

  Depending on what is installed, default order is: torch -> fftw -> numpy

Returns:

spec – 1d Fast-Fourier Transformed data, along axis. Same shape as data, except axis now has length n_fft.

Return type:

ndarray, shape=(…,n_fft,…)

ifft(spec, n_fft, axis=0, fft_method=None)

Compute inverse 1d discrete Fast Fourier transform along given array axis, using given method

Thin wrapper around multiple low-level iFFT implementations, which may be chosen based on performance and availability (in our hands, torch is fastest for virtually all tests).

Parameters:

• spec (ndarray, shape=(...,n_freqs,...)) – Spectral data to compute inverse 1d FFT on. Shape is arbitrary; iFFT is performed along axis (corresponding to frequency), independently across any/all other array axes (which might be data channels, trials, etc.).

• n_fft (int, default: spec.shape[axis]) – Length of the transformed axis of the output. If smaller than the length of the input, the input is cropped. If larger, the input is zero-padded with zeros. If not given, the length of the input along the axis specified by axis is used.

• axis (int, default: 0 (1st array axis)) – Axis over which to compute the FFT. If not given, the first axis is used.

• fft_method (str, default: 'torch' (if available)) –

  Which underlying low-level inverse FFT implementation to use. Options:

  • ’torch’Torch’s FFT on the GPU, using :func:torch.fft.ifft.

    Across a range of tests, this is by far the fastest method.

  • ’fftw’FFTW’s optimized FFT, using :func:pyfftw.interfaces.scipy_fftpack.ifft

    FFTW is faster than scipy/numpy for data length != power of 2 (but slower than torch).

  • ’numpy’ : Numpy’s FFT implementation, using :func:np.fft.ifft

  • ’scipy’ : Scipy’s FFT implementation, using :func:scipy.fft.ifft

  Depending on what is installed, default order is: torch -> fftw -> numpy

Returns:

data – 1d inverse Fast-Fourier Transformed data, along axis. Same shape as spec, except axis (usually corresponding to time) now has length n_fft.

Return type:

ndarray, shape=(…,n_fft,…)

complex_to_spec_type(data, spec_type)

Converts complex spectral data to given spectral signal type

Parameters:

• data (ndarray, shape=Any, dtype=complex) – Complex spectral (or time-frequency) data. Arbitrary shape.

• spec_type ({'power','phase','magnitude','real','imag'}) –

  Type of spectral signal to return:

  • ’power’ Spectral power of data

  • ’phase’ Phase of complex spectral data (in radians)

  • ’magnitude’ Magnitude (square root of power) of complex data = signal envelope

  • ’real’ Real part of complex data

  • ’imag’ Imaginary part of complex data

Returns:

data – Computed spectral signal. Same shape as input.

Return type:

ndarray, shape=Any, dtype=complex

power(data)

Compute power from complex spectral data

magnitude(data)

Compute magnitude (square root of power) from complex spectral data

phase(data)

Compute phase of complex spectral data

real(data)

Return real part of complex spectral data

imag(data)

Return imaginary part of complex spectral data

one_sided_to_two_sided(data, freqs, smp_rate, axis=0)

Convert a one-sided Fourier/wavelet transform output to the two-sided equivalent.

Assumes conjugate symmetry across positive and negative frequencies (as is the case only when the original raw signals were real).

Also extrapolates values for f=0, as is necessary for wavelet transforms.

Parameters:

• data (ndarray, shape=(...,n_freqs,...), dtype=complex) – Complex (1-sided) frequency-transformed data. Any arbitary shape.

• freqs (array-like, shape=(n_freqs,)) – Frequency sampling in data

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

• axis (int, default: 0 (1st axis)) – Data axis corresponding to frequency

Returns:

• data (ndarray, shape=(…,2*n_freqs+1,…), dtype=complex) – 2-sided equivalent of input data

• freqs (ndarray, shape=(2*n_freqs+1,)) – List of (positive and negative) freqs in 2-sided output data

simulate_oscillation(frequency, amplitude=5.0, phase=0, noise=1.0, n_trials=1000, freq_sd=0, amp_sd=0, phase_sd=0, smp_rate=1000, time_range=1.0, burst_rate=0, burst_width=4, seed=None)

Generate synthetic data with oscillation at given parameters.

Generate multiple trials with constant oscillatory signal + random additive Gaussian noise.

Parameters:

• frequency (scalar) – Frequency to simulate oscillation at (Hz)

• amplitude (scalar, default: 5.0) – Amplitude of simulated oscillation (a.u.)

• phase (scalar, default: 0) – Phase of oscillation (rad)

• noise (scalar, default: 1.0) – Amplitude of additive Gaussian noise (a.u)

• n_trials (int, default: 1000) – Number of trials/observations to simulate

• freq_sd (scalar, Default: 0 (no inter-trial variation)) – Inter-trial variation in frequency/amplitude/phase, given as Gaussian SD (same units as base parameters, which are used as Gaussian mean)

• amp_sd (scalar, Default: 0 (no inter-trial variation)) – Inter-trial variation in frequency/amplitude/phase, given as Gaussian SD (same units as base parameters, which are used as Gaussian mean)

• phase_sd (scalar, Default: 0 (no inter-trial variation)) – Inter-trial variation in frequency/amplitude/phase, given as Gaussian SD (same units as base parameters, which are used as Gaussian mean)

• smp_rate (int, default: 1000) – Sampling rate for simulated data (Hz)

• time_range (scalar, default: 1 s) – Full time range to simulate oscillation over (s)

• burst_rate (scalar, default: 0 (not bursty)) – Oscillatory burst rate (bursts/trial). Set=0 to simulate constant, non-bursty oscillation.

• burst_width (scalar, default: 4) – Half-width of oscillatory bursts (Gaussian SD, in cycles)

• seed (int, default: None) – Random generator seed for repeatable results. Set=None for unseeded random numbers.

Returns:

data – Simulated oscillation-in-noise data

Return type:

ndarray, shape=(n_timepts,n_trials)