elephant.spectral.welch_psd

elephant.spectral.welch_psd(signal, n_segments=8, len_segment=None, frequency_resolution=None, overlap=0.5, fs=1.0, window='hann', nfft=None, detrend='constant', return_onesided=True, scaling='density', axis=-1)[source]

Estimates power spectrum density (PSD) of a given neo.AnalogSignal using Welch’s method.

The PSD is obtained through the following steps:

  1. Cut the given data into several overlapping segments. The degree of overlap can be specified by parameter overlap (default is 0.5, i.e. segments are overlapped by the half of their length). The number and the length of the segments are determined according to the parameters n_segments, len_segment or frequency_resolution. By default, the data is cut into 8 segments;

  2. Apply a window function to each segment. Hanning window is used by default. This can be changed by giving a window function or an array as parameter window (see Notes [2]);

  3. Compute the periodogram of each segment;

  4. Average the obtained periodograms to yield PSD estimate.

Parameters:
signalneo.AnalogSignal or pq.Quantity or np.ndarray

Time series data, of which PSD is estimated. When signal is pq.Quantity or np.ndarray, sampling frequency should be given through the keyword argument fs. Otherwise, the default value is used (fs = 1.0).

n_segmentsint, optional

Number of segments. The length of segments is adjusted so that overlapping segments cover the entire stretch of the given data. This parameter is ignored if len_segment or frequency_resolution is given. Default: 8

len_segmentint, optional

Length of segments. This parameter is ignored if frequency_resolution is given. If None, it will be determined from other parameters. Default: None

frequency_resolutionpq.Quantity or float, optional

Desired frequency resolution of the obtained PSD estimate in terms of the interval between adjacent frequency bins. When given as a float, it is taken as frequency in Hz. If None, it will be determined from other parameters. Default: None

overlapfloat, optional

Overlap between segments represented as a float number between 0 (no overlap) and 1 (complete overlap). Default: 0.5 (half-overlapped)

fspq.Quantity or float, optional

Specifies the sampling frequency of the input time series. When the input is given as a neo.AnalogSignal, the sampling frequency is taken from its attribute and this parameter is ignored. Default: 1.0

windowstr or tuple or np.ndarray, optional

Desired window to use. See Notes [2]. Default: ‘hann’

nfftint, optional

Length of the FFT used. See Notes [2]. Default: None

detrendstr or function or False, optional

Specifies how to detrend each segment. See Notes [2]. Default: ‘constant’

return_onesidedbool, optional

If True, return a one-sided spectrum for real data. If False return a two-sided spectrum. See Notes [2]. Default: True

scaling{‘density’, ‘spectrum’}, optional

If ‘density’, computes the power spectral density where Pxx has units of V**2/Hz. If ‘spectrum’, computes the power spectrum where Pxx has units of V**2, if signal is measured in V and fs is measured in Hz. See Notes [2]. Default: ‘density’

axisint, optional

Axis along which the periodogram is computed. See Notes [2]. Default: last axis (-1)

Returns:
freqspq.Quantity or np.ndarray

Frequencies associated with the power estimates in psd. freqs is always a vector irrespective of the shape of the input data in signal. If signal is neo.AnalogSignal or pq.Quantity, a pq.Quantity array is returned. Otherwise, a np.ndarray containing frequency in Hz is returned.

psdpq.Quantity or np.ndarray

PSD estimates of the time series in signal. If signal is neo.AnalogSignal, a pq.Quantity array is returned. Otherwise, the return is a np.ndarray.

Raises:
ValueError

If overlap is not in the interval [0, 1).

If frequency_resolution is not positive.

If frequency_resolution is too high for the given data size.

If frequency_resolution is None and len_segment is not a positive number.

If frequency_resolution is None and len_segment is greater than the length of data at axis.

If both frequency_resolution and len_segment are None and n_segments is not a positive number.

If both frequency_resolution and len_segment are None and n_segments is greater than the length of data at axis.

See also

scipy.signal.welch
welch_cohere

Notes

  1. The computation steps used in this function are implemented in scipy.signal module, and this function is a wrapper which provides a proper set of parameters to scipy.signal.welch function.

  2. The parameters window, nfft, detrend, return_onesided, scaling, and axis are directly passed to the scipy.signal.welch function. See the respective descriptions in the docstring of scipy.signal.welch for usage.

  3. When only n_segments is given, parameter nperseg of scipy.signal.welch function is determined according to the expression

    signal.shape[axis] / (n_segments - overlap * (n_segments - 1))

    converted to integer.

Examples

>>> import neo
>>> import numpy as np
>>> import quantities as pq
>>> from elephant.spectral import welch_psd
>>> signal = neo.AnalogSignal(np.cos(np.linspace(0, 2 * np.pi, num=100)),
...     sampling_rate=20 * pq.Hz, units='mV')

Sampling frequency will be taken as signal.sampling_rate.

>>> freq, psd = welch_psd(signal)
>>> freq
array([ 0.        ,  0.90909091,  1.81818182,  2.72727273,  3.63636364,
        4.54545455,  5.45454545,  6.36363636,  7.27272727,  8.18181818,
        9.09090909, 10.        ]) * Hz
>>> psd # noqa
array([[1.09566410e-03, 2.33607943e-02, 1.35436832e-03, 6.74408723e-05,
        1.00810196e-05, 2.40079315e-06, 7.35821437e-07, 2.58361700e-07,
        9.44183422e-08, 3.14573483e-08, 6.82050475e-09, 1.18183354e-10]]) * mV**2/Hz