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:
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;
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]);
Compute the periodogram of each segment;
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
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.
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.
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