Skip to content

Power Spectrum

openseize.spectra.estimators.psd(data, fs, axis=-1, resolution=0.5, window='hann', overlap=0.5, detrend='constant', scaling='density')

A power spectrum (density) estimator using Welch's method.

Welch's method divides data into overlapping segments and averages the modified periodograms for each segment. This estimator can process ndarrays or producer of ndarrays allowing it to estimate the PSD for large data sets.

Parameters:

Name Type Description Default
data Union[npt.NDArray[np.float64], Producer]

An ndarray or producer whose power spectral is to be estimated.

 required
fs float

The sampling rate of the data.

 required
axis int

The sample axis of the data. The estimate will be carried out along this axis.

 -1
resolution float

The frequency resolution of the estimate in Hz. The resolution determines the number of DFT frequencies between [0, fs) to use to estimate the power spectra. The number of DFT points is fs // resolution. The default resolution is 0.5 Hz.

 0.5
window str

A string name for a scipy window applied to each data segment before computing the periodogram of that segment. For a full list of windows see scipy.signal.windows.

 'hann'
overlap float

A percentage in [0, 1) of the data segment overlap . If 0 this estimate is equivalent to Bartletts method (2).

 0.5
detrend str

A string indicating whether to subtract the mean ('constant') or subtract a linear fit ('linear') from each segment prior to computing the estimate for a segment.

 'constant'
scaling str

Determines the normalization to apply to the estimate. If 'spectrum' the estimate will have units V^2 and if 'density' V^2 / Hz.

 'density'

Returns:

Type Description
Tuple[int, npt.NDArray[np.float64], npt.NDArray[np.float64]]

A tuple (n, frequencies, estimate) where:

  • n is the integer number of windows averaged to estimate the PSD
  • frequencies is a 1D array at which the PSD was estimated
  • estimate is a 2-D array of PSD estimates one per channel.

Examples:

>>> # import demo data and make a producer
>>> from openseize.demos import paths
>>> from openseize.file_io.edf import Reader
>>> from openseize import producer
>>> from openseize.spectra.estimators import psd
>>> import matplotlib.pyplot as plt
>>> fp = paths.locate('recording_001.edf')
>>> reader = Reader(fp)
>>> pro = producer(reader, chunksize=10e4, axis=-1)
>>> # Compute the PSD
>>> n, freqs, estimate = psd(pro, fs=5000, axis=-1)
>>> # plot the channel 0 psd
>>> plt.plot(freqs, estimate[0])

References:

  • (1) P. Welch, "The use of the fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms", IEEE Trans. Audio Electroacoust. vol.15, pp. 70-73, 1967

  • (2) M.S. Bartlett, "Periodogram Analysis and Continuous Spectra", Biometrika, vol. 37, pp. 1-16, 1950.

  • (3) B. Porat, "A Course In Digital Signal Processing" Chapters 4 & 13. Wiley and Sons 1997.

Source code in openseize/spectra/estimators.py 
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
def psd(data: Union[npt.NDArray[np.float64], Producer],
        fs: float,
        axis: int = -1,
        resolution: float = 0.5,
        window: str = 'hann',
        overlap: float = 0.5,
        detrend: str = 'constant',
        scaling: str = 'density'
) -> Tuple[int, npt.NDArray[np.float64], npt.NDArray[np.float64]]:
    """A power spectrum (density) estimator using Welch's method.

    Welch's method divides data into overlapping segments and averages the
    modified periodograms for each segment. This estimator can process
    ndarrays or producer of ndarrays allowing it to estimate the PSD for
    large data sets.

    Args:
        data:
            An ndarray or producer whose power spectral is to be estimated.
        fs:
            The sampling rate of the data.
        axis:
            The sample axis of the data. The estimate will be carried out
            along this axis.
        resolution:
            The frequency resolution of the estimate in Hz. The resolution
            determines the number of DFT frequencies between [0, fs) to use
            to estimate the power spectra. The number of DFT points is
            fs // resolution. The default resolution is 0.5 Hz.
        window:
            A string name for a scipy window applied to each data segment
            before computing the periodogram of that segment. For a full
            list of windows see scipy.signal.windows.
        overlap:
            A percentage in [0, 1) of the data segment overlap . If 0 this
            estimate is equivalent to Bartletts method (2).
        detrend:
            A string indicating whether to subtract the mean ('constant') or
            subtract a linear fit ('linear') from each segment prior to
            computing the estimate for a segment.
        scaling:
            Determines the normalization to apply to the estimate. If
            'spectrum' the estimate will have units V^2 and if 'density'
            V^2 / Hz.

    Returns:
        A tuple (n, frequencies, estimate) where:

         - n is the integer number of windows averaged to estimate the PSD
         - frequencies is a 1D array at which the PSD was estimated
         - estimate is a 2-D array of PSD estimates one per channel.

    Examples:
        >>> # import demo data and make a producer
        >>> from openseize.demos import paths
        >>> from openseize.file_io.edf import Reader
        >>> from openseize import producer
        >>> from openseize.spectra.estimators import psd
        >>> import matplotlib.pyplot as plt
        >>> fp = paths.locate('recording_001.edf')
        >>> reader = Reader(fp)
        >>> pro = producer(reader, chunksize=10e4, axis=-1)
        >>> # Compute the PSD
        >>> n, freqs, estimate = psd(pro, fs=5000, axis=-1)
        >>> # plot the channel 0 psd
        >>> plt.plot(freqs, estimate[0])

    References:

    - (1) P. Welch, "The use of the fast Fourier transform for the
        estimation of power spectra: A method based on time averaging
        over short, modified periodograms", IEEE Trans. Audio
        Electroacoust. vol.15, pp. 70-73, 1967

    - (2) M.S. Bartlett, "Periodogram Analysis and Continuous
        Spectra", Biometrika, vol. 37, pp. 1-16, 1950.

    - (3) B. Porat, "A Course In Digital Signal Processing" Chapters
        4 & 13. Wiley and Sons 1997.
    """

    pro = producer(data, chunksize=int(fs), axis=axis)

    # convert requested resolution to DFT pts
    nfft = int(fs / resolution)

    # build a producer of psd estimates, one per welch segment & store
    freqs, psd_pro = nm.welch(pro, fs, nfft, window, overlap, axis, detrend,
                              scaling)

    # compute the average PSD estimate
    result = 0
    for cnt, arr in enumerate(psd_pro, 1):
        result = result + 1 / cnt * (arr - result)

    # pylint misses the cnt variable here
    #pylint: disable-next=undefined-loop-variable
    return cnt, freqs, result #type: ignore