Module `thunderlab.powerspectrum`

Powerspectra and spectrograms for a given frequency resolution

Computation of nfft

next_power_of_two(): round an integer up to the next power of two.
nfft(): compute nfft based on a given frequency resolution.

Decibel

decibel(): transform power to decibel.
power(): transform decibel to power.

Power spectra

psd(): power spectrum for a given frequency resolution.
multi_psd(): power spectra for consecutive data windows and mutiple frequency resolutions.
spectrogram(): spectrogram of a given frequency resolution and overlap fraction.

Power spectrum analysis

peak_freqs(): peak frequencies computed from power spectra of data snippets.

Visualization

plot_decibel_psd(): plot power spectrum in decibel.

Configuration parameter

add_multi_psd_config(): add parameters for multi_psd() to configuration.
multi_psd_args(): retrieve parameters for mulit_psd() from configuration.

Functions

def next_power_of_two(n)

The next integer power of two.

Parameters

n : int: A positive number.

Returns

m : int: The next integer power of two equal or larger than n.

def nfft(rate, freq_resolution, min_nfft=16, max_nfft=None)

Required number of samples for an FFT of a given frequency resolution.

Note that the returned number of FFT samples results in frequency intervals that are smaller or equal to freq_resolution.

Parameters

rate : float: Sampling rate of the data in Hertz.
freq_resolution : float: Minimum frequency resolution in Hertz.
min_nfft : int: Smallest value of nfft to be used.
max_nfft : int or None: If not None, largest value of nfft to be used.

Returns

nfft : int: Number of FFT points.

def decibel(power, ref_power=1.0, min_power=1e-20)

Transform power to decibel relative to ref_power.

$decibel = 10 \cdot \log_{10}(power/ref\_power)$ Power values smaller than min_power are set to -np.inf.

Parameters

power : float or array: Power values, for example from a power spectrum or spectrogram.
ref_power : float or None or 'peak': Reference power for computing decibel. If set to None or 'peak', the maximum power is used.
min_power : float: Power values smaller than min_power are set to -np.inf.

Returns

decibel_psd : array: Power values in decibel relative to ref_power.

def power(decibel, ref_power=1.0)

Transform decibel back to power relative to ref_power.

$power = ref\_power \cdot 10^{decibel/10}$

Parameters

decibel : array: Decibel values of the power spectrum or spectrogram.
ref_power : float: Reference power for computing power.

Returns

power : array: Power values of the power spectrum or spectrogram.

def psd(data, ratetime, freq_resolution, min_nfft=16, max_nfft=None, overlap_frac=0.5, detrend='constant', window='hann')

Power spectrum density of a given frequency resolution.

NFFT is computed from the requested frequency resolution and the sampling rate. Check the returned frequency array for the actually used frequency resolution. The frequency intervals are smaller or equal to freq_resolution. NFFT can be retrieved by dividing the sampling rate by the actual frequency resolution:

freq, power = psd(data, samplingrate, 0.1)
df = np.mean(np.diff(freq))  # the actual frequency resolution
nfft = int(samplingrate/df)

Uses scipy signal.welch() if available, otherwise matplotlib.mlab.psd().

Parameters

data : 1-D array: Data from which power spectra are computed.
ratetime : float or array: If float, sampling rate of the data in Hertz. If array, assume ratetime to be the time array corresponding to the data. Compute sampling rate as 1/(ratetime[1]-ratetime[0]).
freq_resolution : float: Frequency resolution of the psd in Hertz.
min_nfft : int: Smallest value of nfft to be used.
max_nfft : int or None: If not None, largest value of nfft to be used.
overlap_frac : float: Fraction of overlap for the fft windows.
detrend : string: If 'constant' or 'mean' subtract mean of data. If 'linear' subtract line fitted to the data. If 'none' do not detrend the data.
window : string: Function used for windowing data segements. One of hann, blackman, hamming, bartlett, boxcar, triang, parzen, bohman, blackmanharris, nuttall, fattop, barthann (see scipy.signal window functions).

Returns

freq : 1-D array: Frequencies corresponding to power array.
power : 1-D array: Power spectral density in [data]^2/Hz.

def multi_psd(data, ratetime, freq_resolution=0.2, num_windows=1, min_nfft=16, overlap_frac=0.5, detrend='constant', window='hann')

Power spectra computed for consecutive data windows.

See also psd() for more information on power spectra with given frequency resolution.

Parameters

data : 1-D array: Data from which power spectra are computed.
ratetime : float or array: If float, sampling rate of the data in Hertz. If array, assume ratetime to be the time array corresponding to the data. Compute sampling rate as 1/(ratetime[1]-ratetime[0]).
freq_resolution : float: Frequency resolution of psd in Hertz.
num_windows : int: Data are chopped into num_windows segments that overlap by half for which power spectra are computed.
min_nfft : int: Smallest value of nfft to be used.
overlap_frac : float: Fraction of overlap for the fft windows within a single power spectrum.
detrend : string: If 'constant' subtract mean of data. If 'linear' subtract line fitted to the data. If 'none' do not deternd the data.
window : string: Function used for windowing data segements. One of hann, blackman, hamming, bartlett, boxcar, triang, parzen, bohman, blackmanharris, nuttall, fattop, barthann (see scipy.signal window functions).

Returns

multi_psd_data : list of 2-D arrays: List of the power spectra for each window and frequency resolution (psd_data[i][freq, power()]).

def spectrogram(data, ratetime, freq_resolution=0.2, min_nfft=16, max_nfft=None, overlap_frac=0.5, detrend='constant', window='hann')

Spectrogram of a given frequency resolution.

Check the returned frequency array for the actually used frequency resolution. The actual frequency resolution is smaller or equal to freq_resolution. The used number of data points per FFT segment (NFFT) is the sampling rate divided by the actual frequency resolution:

spec, freq, time = spectrum(data, samplingrate, 0.1) # request 0.1Hz resolution
df = np.mean(np.diff(freq))  # the actual frequency resolution
nfft = int(samplingrate/df)

Parameters

data : 1D or 2D array of floats: Data for the spectrograms. First dimension is time, optional second dimension is channel.
ratetime : float or array: If float, sampling rate of the data in Hertz. If array, assume ratetime to be the time array corresponding to the data. The sampling rate is then computed as 1/(ratetime[1]-ratetime[0]).
freq_resolution : float: Frequency resolution for the spectrogram in Hertz. See nfft() for details.
min_nfft : int: Smallest value of nfft to be used. See nfft() for details.
max_nfft : int or None: If not None, largest value of nfft to be used. See nfft() for details.
overlap_frac : float: Overlap of the nffts (0 = no overlap; 1 = complete overlap).
detrend : string or False: If 'constant' subtract mean of each data segment. If 'linear' subtract line fitted to each data segment. If False do not detrend the data segments.
window : string: Function used for windowing data segements. One of hann, blackman, hamming, bartlett, boxcar, triang, parzen, bohman, blackmanharris, nuttall, fattop, barthann, tukey (see scipy.signal window functions).

Returns

freqs : array: Frequencies of the spectrogram.
time : array: Time of the nfft windows.
spectrum : 2D or 3D array: Power spectral density for each frequency and time. First dimension is frequency and second dimension is time. Optional last dimension is channel.

def plot_decibel_psd(ax, freqs, power, ref_power=1.0, min_power=1e-20, log_freq=False, min_freq=0.0, max_freq=2000.0, ymarg=0.0, **kwargs)

Plot the powerspectum in decibel relative to ref_power.

Parameters

ax:
Axis for plot.
freqs : 1-D array: Frequency array of the power spectrum.
power : 1-D array: Power values of the power spectrum.
ref_power : float: Reference power for computing decibel. If set to None the maximum power is used.
min_power : float: Power values smaller than min_power are set to np.nan.
log_freq : boolean: Logarithmic (True) or linear (False) frequency axis.
min_freq : float: Limits of frequency axis are set to (min_freq, max_freq) if max_freq is greater than zero
max_freq : float: Limits of frequency axis are set to (min_freq, max_freq) and limits of power axis are computed from powers below max_freq if max_freq is greater than zero
ymarg : float: Add this to the maximum decibel power for setting the ylim.
kwargs : dict: Plot parameter that are passed on to the plot() function.

def peak_freqs(onsets, offsets, data, rate, freq_resolution=0.2, thresh=None, **kwargs)

Peak frequencies computed from power spectra of data snippets.

Parameters

onsets : array of ints: Indices indicating the onsets of the snippets in data.
offsets : array of ints: Indices indicating the offsets of the snippets in data.
data : 1-D array: Data array that contains the data snippets defined by onsets and offsets.
rate : float: Sampling rate of data in Hertz.
freq_resolution : float: Desired frequency resolution of the computed power spectra in Hertz.
thresh : None or float: If not None than this is the threshold required for the minimum height of the peak in the decibel power spectrum. If the peak is too small than the peak frequency of that snippet is set to NaN.
kwargs : dict: Further arguments passed on to psd().

Returns

freqs : array of floats: For each data snippet the frequency of the maximum power.

def add_multi_psd_config(cfg, freq_resolution=0.2, num_resolutions=1, num_windows=1)

Add all parameters needed for the multi_psd() function as a new section to a configuration.

Parameters

cfg : ConfigFile: The configuration.

See multi_psd() for details on the remaining arguments.

def multi_psd_args(cfg)

Translates a configuration to the respective parameter names of the multi_psd() function.

The return value can then be passed as key-word arguments to this functions.

Parameters

cfg : ConfigFile: The configuration.

Returns

a : dict: Dictionary with names of arguments of the multi_psd() function and their values as supplied by cfg.

def main()