psiaudio.util module

analyze_mic_sens(ref_waveforms, exp_waveforms, vrms, ref_mic_gain, exp_mic_gain, output_gain, ref_mic_sens, **kwargs)

analyze_tone(waveforms, frequency, fs, mic_gain, trim=0, thd_harmonics=3)

as_numeric(x)

Ensure input is a numeric numpy array or pandas block.

Parameters:

x (array_like) – Input data.

Returns:

Numeric array data.

Return type:

array_like

band_to_spectrum_level(band_db, n)

Convert overall band level to spectrum level

Parameters:

• band_db (float) – Band level in dB re. reference (e.g., dB SPL).

• n (int) – Number of bands. If you are trying to compare this to an FFT with 1 Hz resolution, this is equal to the bandwidth (in Hz). If you are comparing this to an FFT with a bin size of 500 Hz, then this is equal to the bandwidth / 500.

Returns:

spectrum_db – Spectrum level in same units as band level.

Return type:

float

bin_array(number, bits)

Return binary representation of an integer as an integer array >>> bin_array(8, 4) [0, 0, 0, 1] >>> bin_array(3, 4) [1, 1, 0, 0]

NOTE: This function has not been profiled for speed.

check_interleaved_octaves(freqs, min_octaves=1)

Ensure that frequencies are spaced at least an octave apart

Parameters:

• freqs (ordered sequence) – Sequence of frequencies in the desired ordering.

• min_octaves (float) – Minimum octave spacing to enforce (this is multiplied by 0.95 to allow for some fudge factor if frequencies were rounded).

Note

• If you are rounding frequencies to the nearest Hz, the actual octaves spacing may be slightly less or more than the desired spacing due to the rounding. We multiply min_octaves by 0.99 to allow for this small error in octave spacing.

• This also checks that the last frequency is an octave from the first frequency.

csd(s, window=None, detrend='linear')

Compute cross-spectral density using rfft.

Parameters:

• s (array_like) – Input signal.

• window (str or None, optional) – Window to apply to the signal, by default None.

• detrend (str or None, optional) – Detrending type, by default ‘linear’.

Returns:

Cross-spectral density.

Return type:

ndarray

csd_df(s, fs, *args, **kw)

csd_to_signal(csd)

Convert cross-spectral density back to signal.

Parameters:

csd (ndarray) – Cross-spectral density.

Returns:

Time-domain signal.

Return type:

ndarray

db(target, reference=1)

Convert target amplitude to decibels relative to a reference.

Parameters:

• target (array_like) – Target amplitude.

• reference (float or array_like, optional) – Reference amplitude, by default 1.

Returns:

Decibels (dB) computed as 20 * log10(target / reference).

Return type:

array_like

dbi(db, reference=1)

Inverse decibel conversion.

Parameters:

• db (array_like) – Array of decibel values.

• reference (float or array_like, optional) – Reference amplitude, by default 1.

Returns:

Linear amplitude.

Return type:

array_like

dbtopa(db)

Convert dB SPL to Pascal

\[10^{dB/20.0}/(20\cdot10^{-6})\]

>>> round(dbtopa(94), 4)
1.0024
>>> dbtopa(100)
2.0
>>> dbtopa(120)
20.0
>>> patodb(dbtopa(94.0))
94.0

Will also take sequences:

>>> print(dbtopa([80, 100, 120]))
[ 0.2  2.  20. ]

debounce_epochs(epochs, debounce)

Given a 2D array of epochs in the format [[start time, end time], …], throw out all epochs that are shorter than the minimum sample duration. After discarding these epochs, combine remaining epochs if they are within the minimum sample duration of each other.

diff_matrix(n_chan, reference, labels=None)

edge_falling(TTL)

edge_rising(TTL)

epochs(x, pad=0)

Given a boolean array, where 1 = epoch, return indices of epochs (first column is the index where x goes from 0 to 1 and second column is index where x goes from 1 to 0.

epochs_contain(epochs, ts)

Returns True if ts falls within one of the epoch boundaries Epochs must be sorted.

epochs_overlap(a, b)

Returns True where b falls within boundaries of epoch in a Epochs must be sorted.

freq_smooth(frequency, power, bandwidth=20)

Uses Konno & Ohmachi (1998) algorithm

get_cb(cb, suffix=None)

Create a callback function to indicate progress.

Parameters:

• cb ({'tqdm', None}) – If ‘tqdm’, creates a text-based progressbar. If None, no progress is reported.

• suffix (str, optional) – Suffix to append to the progress bar message.

Returns:

Function that can be called iteratively with a fractional completion value.

Return type:

callable

golay_ir(n, a, b, a_signal, b_signal)

Estimate system impulse response from Golay sequence

Implements algorithm described in Zhou et al. 1992

golay_pair(n=15)

Generate pair of Golay sequences

golay_tf(a, b, a_signal, b_signal, fs)

Estimate system transfer function from Golay sequence

Implements algorithm as described in Zhou et al. 1992.

iir(psd, phase, frequency, cutoff=None, phase_correction=None, truncate=None, truncate_spectrum=False, reference='mean')

Given the impulse response, compute the inverse impulse response.

Parameters:

TODO (#)

Note

Specifying the cutoff range is highly recommended to get a well-behaved function.

int_to_TTL(a, width)

Converts a 1D array of integers to a 2D boolean array based on the binary representation of each integer. Primarily used in conjunction with TDT’s FromBits component to reduce the overhead of storing and transferring TTL data. Since a TTL can be represented as a single bit (0 or 1), it is wasteful to cast the TTL to an int32 before storing the data in a buffer. FromBits combines up to 6 TTL channels into a single int32 word. Since only the first 6 bits are used to store the 6 TTL channels, the data can be reduced further:

1. Using ShufTo8, 24 TTL channels can be stored in a single index of a serial buffer.

2. Using CompTo8, 4 consecutive samples of data from 6 TTL channels can be stored in a single index of a serial buffer.

3. Combining ShufTo16 and CompTo16, store 2 consecutive samples of data from 12 TTL channels in a single index of a serial buffer. Using this approach, the memory overhead and amount of data being transferred has been reduced by a factor of 24. This function uses Numpy’s bitshift and bitmask operators, so the algorithm should be pretty efficient.

Parameters:

a (array_like) – Sequence of integers to expand into the corresponding boolean array. The dtype (either int8, int16 or int32) of the array is used to figure out the size of the second dimension. This will depend on your combination of FromBits and the shuffle/compression components.

Returns:

bitfield – 2D boolean array repesenting the bits in little-endian order

Return type:

array

Example

>>> int_to_TTL([4, 8, 5], width=6).T
array([[False, False,  True, False, False, False],
       [False, False, False,  True, False, False],
       [ True, False,  True, False, False, False]])

interleave_octaves(freqs, min_octaves=1)

Return correct ordering for frequencies in interleaved paradigm as per. Buran et al.

This function works with both kHz and Hz.

>>> interleave_octaves([2, 2.8, 4, 5.6, 8])
[8, 4, 2, 5.6, 2.8]
>>> interleave_octaves([2000, 2800, 4000, 5600, 8000])
[8000, 4000, 2000, 5600, 2800]

If a set of frequencies cannot appropriately be ordered, a ValueError is raised. In this example, the first and last frequences are within one octave.

>>> interleave_octaves([2000, 2800, 4000])
Traceback (most recent call last):
  ...
ValueError: Unable to interleave 4000, 2000, 2800 so that adjacent
frequencies are spaced at least 1.0 octaves apart. This can often be fixed
by increasing the range of frequencies tested. For example, if you are
assessing frequencies spaced 0.5 octaves, you need at least four
frequencies to be able order them so that adjacent frequencies are at least
1.0 octaves apart.

You can use a different octave spacing.

>>> interleave_octaves([2000, 2800, 4000], 0.5)
[4000, 2800, 2000]

ir_iir(impulse_response, fs, smooth=None, *args, **kwargs)

nearest_octave(x, step, si_prefix='')

Coerces x to the nearest octave step with some sensible defaults

Parameters:

• x ({int, float}) – Value to coerce.

• step (float) – Octave step size to coerce to.

• {'' (si_prefix) – If ‘’, x is assumed to be in Hz. If ‘k’, x is assumed to be in kHz. This is important because the standard coercion algorithm will coerce 8000 Hz to 8192 Hz otherwise.

• 'k'} – If ‘’, x is assumed to be in Hz. If ‘k’, x is assumed to be in kHz. This is important because the standard coercion algorithm will coerce 8000 Hz to 8192 Hz otherwise.

Examples

>>> nearest_octave(45200, 0.5)
45255
>>> nearest_octave(11000, 0.5)
11314
>>> nearest_octave(8000, 0.5)
8000
>>> nearest_octave(45.2, 0.5, 'k')
45.255
>>> nearest_octave(11, 0.5, 'k')
11.314
>>> nearest_octave(8.0, 0.5, 'k')
8.0

normalize_rms(waveform, out=None)

Normalize RMS power to 1.

Typically used when generating a noise waveform that will be scaled by a calibration factor.

Parameters:

• waveform (array_like) – Input array.

• out (array_like, optional) – An array to store the output. Must be the same shape as waveform.

Returns:

RMS normalized waveform.

Return type:

array_like

octave_band_freqs(fc, octaves)

Calculate frequency limits of a band given center frequency and width in octaves.

The center frequency is the geometric mean of the upper and lower frequency limits.

Parameters:

• fc (float) – Center frequency of band.

• octaves (float) – Bandwidth in octaves.

Returns:

tuple

Return type:

A tuple containing the lower and upper frequencies.

Examples

>>> fl, fh = octave_band_freqs(8e3, 1/8)
>>> print(f'{fl:.0f} to {fh:.0f}'
7660 to 8354

>>> fl, fh = octave_band_freqs(8, 1/8)
>>> print(f'{fl:.2f} to {fh:.2f}'
7.66 to 8.35

octave_space(lb, ub, step, mode='nearest')

>>> print(octave_space(4, 32, 1.0))
[ 4.  8. 16. 32.]

>>> freq = octave_space(0.5, 50.0, 0.25, 'nearest')
>>> print(round(min(freq), 2))
0.5
>>> print(round(max(freq), 2))
53.82

>>> freq = octave_space(0.5, 50.0, 0.25, 'bounded')
>>> print(round(min(freq), 2))
0.5
>>> print(round(max(freq), 2))
45.25

patodb(pa)

Convert Pascal to dB SPL

\[20*log10(pa/20e-6)\]

>>> round(patodb(1))
94
>>> patodb(2)
100.0
>>> patodb(0.2)
80.0

Will also take sequences: >>> print(patodb([0.2, 2.0, 20.0])) [ 80. 100. 120.]

phase(s, fs, window=None, waveform_averages=None, unwrap=True)

Extract phase from a signal.

Parameters:

• s (array_like) – Input signal.

• fs (float) – Sampling rate.

• window (str, optional) – Window name, by default None.

• waveform_averages (int, optional) – Number of averages, by default None.

• unwrap (bool, optional) – Whether to unwrap phase angle, by default True.

Returns:

Phase of the signal in radians.

Return type:

ndarray

phase_df(s, fs, *args, waveform_averages=None, **kw)

process_tone(fs, signal, frequency, min_snr=None, max_thd=None, thd_harmonics=3, silence=None)

Compute the RMS at the specified frequency. Check for distortion.

Parameters:

• fs (float) – Sampling frequency of signal

• signal (ndarray) – Last dimension must be time. If more than one dimension, first dimension must be repetition.

• frequency (float) – Frequency (Hz) to analyze

• min_snr ({None, float}) – If specified, must provide a noise floor measure (silence). The ratio, in dB, of signal RMS to silence RMS must be greater than min_snr. If not, a CalibrationNFError is raised.

• max_thd ({None, float}) – If specified, ensures that the total harmonic distortion, as a percentage, is less than max_thd. If not, a CalibrationTHDError is raised.

• thd_harmonics (int) – Number of harmonics to compute. If you pick too many, some harmonics may be above the Nyquist frequency and you’ll get an exception.

• thd_harmonics – Number of harmonics to compute. If you pick too many, some harmonics may be above the Nyquist frequency and you’ll get an exception.

• silence ({None, ndarray}) – Noise floor measurement. Required for min_snr. Shape must match signal in all dimensions except the first and last.

Returns:

result – Series will be indexed with RMS, SNR, THD and frequency. DataFrame will contain columns for RMS, SNR, THD and frequency. The return type will depend on the dimensionality of the input array.

Return type:

pandas Series or DataFrame

psd(s, fs, window=None, waveform_averages=None, trim_samples=True, **kw)

psd_bootstrap_loop(arrays, fs, n_draw=100, n_bootstrap=100, rng=None, window=None, callback='tqdm', calculate=None)

Calculate the normalized PSD across trials using a boostrapping algorithm

To estmate the noise floor, the CSD of each trial is computed and then the phases of the CSD are randomized.

Parameters:

• arrays (array or list of arrays) – Signal to compute PSD noise floor over. Must be trial x time. If a list of arrays is provided (e.g., one for each stimulus polarity), n_draw will be balanced across the arrays.

• fs (float) – Sampling rate of signal

• n_draw (int) – Number of trials to draw on each bootstrap cycle. If None, draw all samples (with replacement). Must be a multiple of the number of arrays.

• n_bootstrap (int) – Number of bootstrap cycles.

• rng (instance of RandomState) – If provided, this will be used to drive the bootstrapping algorithm. Useful when you need to “freeze” results (e.g., for publication).

• window ({None, string}) – Type of window to use when calculating bootstrapped PSD.

• Result

• ------

• psd_bs (DataFrame) – Pandas DataFrame indexed by frequency. Columns include psd_nf, the noise floor as estimated by the bootstrapping algorithm.

Notes

This algorithm is adapted from Zhu et al. 2013 (JASA).

psd_bootstrap_vec(x, fs, n_draw=400, n_bootstrap=100, rng=None, window=None)

Calculate the normalized PSD across trials using a boostrapping algorithm

To estmate the noise floor, the CSD of each trial is computed and then the phases of the CSD are randomized.

Parameters:

• x (array) – Signal to compute PSD noise floor over. Must be trial x time.

• fs (float) – Sampling rate of signal

• n_draw (int) – Number of trials to draw on each bootstrap cycle.

• n_bootstrap (int) – Number of bootstrap cycles.

• rng (instance of RandomState) – If provided, this will be used to drive the bootstrapping algorithm. Useful when you need to “freeze” results (e.g., for publication).

• window ({None, string}) – Type of window to use when calculating bootstrapped PSD.

• Result

• ------

• psd_bs (DataFrame) – Pandas DataFrame indexed by frequency. Columns include psd_nf, the noise floor as estimated by the bootstrapping algorithm.

Notes

This algorithm is adapted from Zhu et al. 2013 (JASA).

psd_df(s, fs, *args, waveform_averages=None, **kw)

Compute PSD and return as a dataframe with columns indexed by frequency

Parameters:

• s ({array, DataFrame}) – Data to compute PSD over. The last axis (columns for DataFrame) is always assumed to be time. If array, only 2D arrays are supported (for conversion into DataFrame).

• fs (float) – Sampling rate of data.

• window ({None, string}) – Applies given window to signal before computing FFT. If None, no windowing is performed. Name of window must be a valid window name that can be passed to scipy.signal.get_window.

• waveform_averages ({None, int}) – Number of chunks to segment time data into before computing PSD. Averaging is done after computing the PSD.

• trim_samples (bool) – If true, remove excess time samples so that waveform can be split into waveform_averages segments of equal size. If False and the number of timepoints is not an integer multiple of waveform_averages, an error will be raised.

:param Additional arguments are passed to psd().:

Returns:

psd – DataFrame with frequency as columns. If s was a DataFrame, the original index of s will be preserved. If s was an array, psd will have a simple integer index.

Return type:

pd.DataFrame

psd_freq(s, fs)

resample_fft(waveform, fs, target_fs)

resample_poly(waveform, fs, target_fs)

rms(s, detrend=False, axis=-1)

rms_rfft(x)

rms_rfft_db(x, *args, **kw)

smooth_epochs(epochs)

Given a 2D array of epochs in the format [[start time, end time], …], identify and remove all overlapping epochs such that:

[ epoch ] [ epoch ]

[ epoch ]

Will become:

[ epoch ] [ epoch ]

Epochs do not need to be ordered when provided; however, they will be returned ordered.

spectrum_to_band_level(spectrum_db, n)

Convert spectrum level to overall level

Parameters:

• spectrum_db (float) – Spectrum level in dB re. reference (e.g., dB SPL).

• n (int) – Number of bands. If you are trying to compare this to an FFT with 1 Hz resolution, this is equal to the bandwidth (in Hz). If you are comparing this to an FFT with a bin size of 500 Hz, then this is equal to the bandwidth / 500.

Returns:

band_db – Band level in same units as spectrum level.

Return type:

float

Examples

Calculate band level for broadband noise of 20 msec duration generated at 100 kHz with a spectrum level of 94 dB SPL. >>> bin_width = 1 / 20e-3 >>> f_max = 100e3 / 2 >>> n_bins = f_max / bin_width >>> spectrum_to_band_level(94, n_bins) 124.0

summarize_golay(fs, a, b, a_response, b_response, waveform_averages=None)

thd(s, fs, frequency, harmonics=3, window=None)

tone_conv(s, fs, frequency, window=None, detrend='linear')

tone_phase_conv(s, fs, frequency, window=None)

tone_phase_fft(s, fs, frequency, window=None)

tone_power_conv(s, fs, frequency, window=None, detrend='linear')

tone_power_conv_nf(s, fs, frequency, window=None)

tone_power_fft(s, fs, frequency, window=None)

transfer_function(stimulus, response, fs)

truncated_ifft(spectrum, original_fs, truncated_fs)

ts(TTL)