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()