pyfar.dsp#

Digital signal processing functions.

Classes:

InterpolateSpectrum(data, method, kind[, ...])

Interpolate an incomplete spectrum to a complete single sided spectrum.

Functions:

`average`(signal[, mode, caxis, weights, ...])	Average multi-channel signals.
`convolve`(signal1, signal2[, mode, method])	Convolve two signals.
`decibel`(signal[, domain, log_prefix, ...])	Convert data of the selected signal domain into decibels (dB).
`deconvolve`(system_output, system_input[, ...])	Calculate transfer functions by spectral deconvolution of two signals.
`energy`(signal)	Computes the channel wise energy in the time domain.
`find_impulse_response_delay`(impulse_response)	Find the delay in sub-sample values of an impulse response.
`find_impulse_response_start`(impulse_response)	Find the start sample of an impulse response.
`fractional_time_shift`(signal, shift[, unit, ...])	Apply fractional time shift to input data.
`group_delay`(signal[, frequencies, method])	Returns the group delay of a signal in samples.
`kaiser_window_beta`(A)	Return a shape parameter beta to create kaiser window based on desired side lobe suppression in dB.
`linear_phase`(signal, group_delay[, unit])	Set the phase to a linear phase with a specified group delay.
`minimum_phase`(signal[, n_fft, truncate])	Calculate the minimum phase equivalent of a finite impulse response.
`normalize`(signal[, reference_method, ...])	Apply a normalization.
`pad_zeros`(signal, pad_width[, mode])	Pad a signal with zeros in the time domain.
`phase`(signal[, deg, unwrap])	Returns the phase for a given signal object.
`power`(signal)	Compute the power of a signal.
`regularized_spectrum_inversion`(signal, ...)	Invert the spectrum of a signal applying frequency dependent regularization.
`resample`(signal, sampling_rate[, ...])	Resample signal to new sampling rate.
`rms`(signal)	Compute the root mean square (RMS) of a signal.
`smooth_fractional_octave`(signal, num_fractions)	Smooth spectrum with a fractional octave width.
`soft_limit_spectrum`(signal, limit, knee[, ...])	Soft limiting the magniude spectrum.
`spectrogram`(signal[, window, window_length, ...])	Compute the magnitude spectrum versus time.
`time_shift`(signal, shift[, mode, unit, ...])	Apply a cyclic or linear time-shift to a signal.
`time_window`(signal, interval[, window, ...])	Apply time window to signal.
`wrap_to_2pi`(x)	Wraps phase to 2 pi.
`zero_phase`(signal)	Calculate zero phase signal.

class pyfar.dsp.InterpolateSpectrum(data, method, kind, fscale='linear', clip=False)[source]#

Bases: object

Interpolate an incomplete spectrum to a complete single sided spectrum.

This is intended to interpolate transfer functions, for example sparse spectra that are defined only at octave frequencies or incomplete spectra from numerical simulations.

Parameters:

data (FrequencyData) – Input data to be interpolated.
method (string) –
Specifies the input data for the interpolation

'complex'
Separate interpolation of the real and imaginary part

'magnitude_phase'
Separate interpolation if the magnitude and unwrapped phase values

'magnitude'
Interpolate the magnitude values only. Results in a zero phase signal, which is symmetric around the first sample. This phase response might not be ideal for many applications. Minimum and linear phase responses can be generated with minimum_phase and linear_phase.
kind (tuple) –
Three element tuple ('first', 'second', 'third') that specifies the kind of inter/extrapolation below the lowest frequency (first), between the lowest and highest frequency (second), and above the highest frequency (third).

The individual strings have to be

'zero', slinear, 'quadratic', 'cubic'
Spline interpolation of zeroth, first, second or third order

'previous', 'next'
Simply return the previous or next value of the point

'nearest-up', 'nearest'
Differ when interpolating half-integers (e.g. 0.5, 1.5) in that 'nearest-up' rounds up and 'nearest' rounds down.

The interpolation is done using scipy.interpolate.interp1d.
fscale (string, optional) –

'linear'
Interpolate on a linear frequency axis.

'log'
Interpolate on a logarithmic frequency axis.

The default is 'linear'.
clip (bool, tuple) – The interpolated magnitude response is clipped to the range specified by this two element tuple. E.g., clip=(0, 1) will assure that no values smaller than 0 and larger than 1 occur in the interpolated magnitude response. The clipping is applied after the interpolation. The default is False which does not clip the data.

Returns:

interpolator – The interpolator can be called to interpolate the data (see examples below). It returns a Signal and has the following parameters

n_samplesint: Length of the interpolated time signal in samples
sampling_rate: int: Sampling rate of the output signal in Hz
showbool, optional: Show a plot of the input and output data. The default is False.

Return type:

InterpolateSpectrum

Examples

Interpolate a magnitude spectrum, add an artificial linear phase and inspect the results. Note that a similar plot can be created by the interpolator object by signal = interpolator(64, 44100, show=True)

>>> import pyfar as pf
>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>>
>>> pf.plot.use()
>>> _, ax = plt.subplots(2, 3)
>>>
>>> # generate data
>>> data = pf.FrequencyData([1, 0], [5e3, 20e3])
>>>
>>> # interpolate and plot
>>> for ff, fscale in enumerate(["linear", "log"]):
>>>     interpolator = pf.dsp.InterpolateSpectrum(
...         data, 'magnitude', ('nearest', 'linear', 'nearest'),
...         fscale)
>>>
>>>     # interpolate to 64 samples linear phase impulse response
>>>     signal = interpolator(64, 44100)
>>>     signal = pf.dsp.linear_phase(signal, 32)
>>>
>>>     # time signal (linear and logarithmic amplitude)
>>>     pf.plot.time(signal, ax=ax[ff, 0], unit='ms', dB=True)
>>>     # frequency plot (linear x-axis)
>>>     pf.plot.freq(
...         signal, dB=False, freq_scale="linear", ax=ax[ff, 1])
>>>     pf.plot.freq(data, dB=False, freq_scale="linear",
...                  ax=ax[ff, 1], c='r', ls='', marker='.')
>>>     ax[ff, 1].set_xlim(0, signal.sampling_rate/2)
>>>     ax[ff, 1].set_title(
...         f"Interpolated on {fscale} frequency scale")
>>>     # frequency plot (log x-axis)
>>>     pf.plot.freq(signal, dB=False, ax=ax[ff, 2], label='input')
>>>     pf.plot.freq(data, dB=False, ax=ax[ff, 2],
...                  c='r', ls='', marker='.', label='output')
>>>     ax[ff, 2].set_xlim(2e3, signal.sampling_rate/2)
>>>     ax[ff, 2].legend(loc='best')

(Source code, png, hires.png, pdf)

pyfar.dsp.average(signal, mode='linear', caxis=None, weights=None, keepdims=False, nan_policy='raise')[source]#

Average multi-channel signals.

Parameters:

signal (Signal, TimeData, FrequencyData) – Input signal.
mode (string) –

'linear'
Average signal.time if the signal is in the time domain and signal.freq if the signal is in the frequency domain. Note that these operations are equivalent for Signal objects due to the linearity of the averaging and the FFT.

'magnitude_zerophase'
Average the magnitude spectra and discard the phase.

'magnitude_phase'
Average the magnitude spectra and the unwrapped phase separatly.

'power'
Average the power spectra \(|X|^2\) and discard the phase. The squaring of the spectra is reversed before returning the averaged signal.

'log_magnitude_zerophase'
Average the logarithmic magnitude spectra using decibel and discard the phase. The logarithm is reversed before returning the averaged signal.

The default is 'linear'
caxis (None, int, or tuple of ints, optional) – Channel axes (caxis) along which the averaging is done. The caxis denotes an axis of the data inside an audio object but ignores the last axis that contains the time samples or frequency bins. The default None averages across all channels.
weights (array like) – Array with channel weights for averaging the data. Must be broadcastable to cshape. The default is None, which applies equal weights to all channels.
keepdims (bool, optional) – If this is True, the axes which are reduced during the averaging are kept as a dimension with size one. Otherwise, singular dimensions will be squeezed after averaging. The default is False.
nan_policy (string, optional) –
Define how to handle NaNs in input signal.

'propagate'
If the input signal includes NaNs, the corresponding averaged output signal value will be NaN.

'omit'
NaNs will be omitted while averaging. For each NaN value, the number of values used for the average operation is also reduced by one. If a signal contains only NaN values in a specific dimensions, the output will be zero. For example if the second sample of a multi channel signal is always NaN, the average will be zero at the second sample.

'raise'
A 'ValueError' will be raised, if the input signal includes NaNs.

The default is 'raise'.

Returns:

averaged_signal – Averaged input Signal.

Return type:

Signal, TimeData, FrequencyData

Notes

The functions linear_phase and minimum_phase can be used to obtain a phase for magnitude spectra after using a mode that discards the phase.

pyfar.dsp.convolve(signal1, signal2, mode='full', method='overlap_add')[source]#

Convolve two signals.

Parameters:

signal1 (Signal) – The first signal
signal2 (Signal) – The second signal. The channel shape (cshape) of this signal must be broadcastable to the cshape of the first signal. The cshape gives the shape of the data inside an audio object but ignores the number of samples or frequency bins.
mode (string, optional) –
A string indicating the size of the output:

'full'
Compute the full discrete linear convolution of the input signals. The output has the length 'signal1.n_samples + signal2.n_samples - 1' (Default).

'cut'
Compute the complete convolution with full and truncate the result to the length of the longer signal.

'cyclic'
The output is the cyclic convolution of the signals, where the shorter signal is zero-padded to fit the length of the longer one. This is done by computing the complete convolution with 'full', adding the tail (i.e., the part that is truncated for mode='cut' to the beginning of the result) and truncating the result to the length of the longer signal.
method (str {'overlap_add', 'fft'}, optional) –
A string indicating which method to use to calculate the convolution:

'overlap_add'
Convolve using the overlap-add algorithm based on scipy.signal.oaconvolve.. (Default)

'fft'
Convolve using FFT based on scipy.signal.fftconvolve.

See Notes for more details.

Returns:

signal – The result of the convolution. The channel dimension (cdim) matches the bigger cdim of the two input signals. The channel dimension gives the number of dimensions of the audio data excluding the last dimension, which is n_samples for time domain objects and n_bins for frequency domain objects.

Return type:

Signal

Notes

The overlap-add method is generally much faster than fft convolution when one signal is much larger than the other, but can be slower when only a few output values are needed or when the signals have a very similar length. For method='overlap_add', integer data will be cast to float.

Examples

Illustrate the different modes.

>>> import pyfar as pf
>>> s1 = pf.Signal([1, 0.5, 0.5], 1000)
>>> s2 = pf.Signal([1,-1], 1000)
>>> full = pf.dsp.convolve(s1, s2, mode='full')
>>> cut = pf.dsp.convolve(s1, s2, mode='cut')
>>> cyc = pf.dsp.convolve(s1, s2, mode='cyclic')
>>> # Plot input and output
>>> with pf.plot.context():
>>>     fig, ax = plt.subplots(2, 1, sharex=True)
>>>     pf.plot.time(s1, ax=ax[0], label='Signal 1', marker='o',
...                  unit='samples')
>>>     pf.plot.time(s2, ax=ax[0], label='Signal 2', marker='o',
...                  unit='samples')
>>>     ax[0].set_title('Input Signals')
>>>     ax[0].legend()
>>>     pf.plot.time(full, ax=ax[1], label='full', marker='o',
...                  unit='samples')
>>>     pf.plot.time(cut, ax=ax[1], label='cut', ls='--',  marker='o',
...                  unit='samples')
>>>     pf.plot.time(cyc, ax=ax[1], label='cyclic', ls=':', marker='o',
...                  unit='samples')
>>>     ax[1].set_title('Convolution Result')
>>>     ax[1].set_ylim(-1.1, 1.1)
>>>     ax[1].legend()

(Source code, png, hires.png, pdf)

pyfar.dsp.decibel(signal, domain='freq', log_prefix=None, log_reference=1, return_prefix=False)[source]#

Convert data of the selected signal domain into decibels (dB).

The converted data is calculated by the base 10 logarithmic scale: data_in_dB = log_prefix * numpy.log10(data/log_reference). By using a logarithmic scale, the deciBel is able to compare quantities that may have vast ratios between them. As an example, the sound pressure in dB can be calculated as followed:

\[L_p = 20\log_{10}\biggl(\frac{p}{p_0}\biggr),\]

where \(20\) is the logarithmic prefix for sound field quantities and \(p_0\) would be the reference for the sound pressure level. A list of commonly used reference values can be found in the ‘log_reference’ parameters section.

Parameters:

signal (Signal, TimeData, FrequencyData) – The signal which is converted into decibel
domain (str) –
The domain, that is converted to decibels:

'freq'
Convert normalized frequency domain data. Signal must be of type ‘Signal’ or ‘FrequencyData’.

'time'
Convert time domain data. Signal must be of type ‘Signal’ or ‘TimeData’.

'freq_raw'
Convert frequency domain data without normalization. Signal must be of type ‘Signal’.

The default is 'freq'.
log_prefix (int) – The prefix for the dB calculation. The default None, uses 10 for signals with 'psd' and 'power' FFT normalization and 20 otherwise.

log_reference (int or float) –

Reference for the logarithm calculation. List of commonly used values:

log_reference	value
Digital signals (dBFs)	1
Sound pressure \(L_p\) (dB)	2e-5 Pa
Voltage \(L_V\) (dBu)	0.7746 volt
Sound intensity \(L_I\) (dB)	1e-12 W/m²
Voltage \(L_V\) (dBV)	1 volt
Electric power \(L_P\) (dB)	1 watt

The default is 1.

return_prefix (bool, optional) – If return_prefix is True, the function will also return the log_prefix value. This can be used to delogrithmize the data. The default is False.

Returns:

decibel (numpy.ndarray) – The given signal in decibel in chosen domain.
log_prefix (int or float) – Will be returned if return_prefix is set to True.

Examples

>>> import pyfar as pf
>>> signal = pf.signals.noise(41000, rms=[1, 1])
>>> decibel_data = decibel(signal, domain='time')

pyfar.dsp.deconvolve(system_output, system_input, fft_length=None, frequency_range=None, **kwargs)[source]#

Calculate transfer functions by spectral deconvolution of two signals.

The transfer function \(H(\omega)\) is calculated by spectral deconvolution (spectral division).

\[H(\omega) = \frac{Y(\omega)}{X(\omega)},\]

where \(X(\omega)\) is the system input signal and \(Y(\omega)\) the system output. Regularized inversion is used to avoid numerical issues in calculating \(X(\omega)^{-1} = 1/X(\omega)\) for small values of \(X(\omega)\) (see regularized_spectrum_inversion). The system response (transfer function) is thus calculated as

\[H(\omega) = Y(\omega)X(\omega)^{-1}.\]

For more information, refer to [1].

Parameters:

system_output (Signal) – The system output signal (e.g., recorded after passing a device under test). The system output signal is zero padded, if it is shorter than the system input signal.
system_input (Signal) – The system input signal (e.g., used to perform a measurement). The system input signal is zero padded, if it is shorter than the system output signal.
frequency_range (tuple, array_like, double) – The upper and lower frequency limits outside of which the regularization factor is to be applied. The default None bypasses the regularization, which might cause numerical instabilities in case of band-limited system_input. Also see regularized_spectrum_inversion.
fft_length (int or None) – The length the signals system_output and system_input are zero padded to before deconvolving. The default is None. In this case only the shorter signal is padded to the length of the longer signal, no padding is applied when both signals have the same length.
kwargs (key value arguments) – Key value arguments to control the inversion of \(H(\omega)\) are passed to to regularized_spectrum_inversion.
freq_range (tuple, array_like, double) – The upper and lower frequency limits outside of which the regularization factor is to be applied. The default None bypasses the regularization, which might cause numerical instabilities in case of band-limited system_input. Also see regularized_spectrum_inversion. 'freq_range' parameter will be deprecated in pyfar 0.8.0 in favor of 'frequency_range'.

Returns:

system_response – The resulting signal after deconvolution, representing the system response (the transfer function). The fft_norm of is set to 'none'.

Return type:

Signal

References

pyfar.dsp.energy(signal)[source]#

Computes the channel wise energy in the time domain.

\[\sum_{n=0}^{N-1}|x[n]|^2=\frac{1}{N}\sum_{k=0}^{N-1}|X[k]|^2,\]

which is equivalent to the frequency domain computation according to Parseval’s theorem [2].

Parameters:: signal (Signal) – The signal to compute the energy from.
Returns:: data – The channel-wise energy of the input signal.
Return type:: numpy.ndarray

Notes

Due to the calculation based on the time data, the returned energy is independent of the signal’s fft_norm. power and rms can be used to compute the power and the rms of a signal.

References

pyfar.dsp.find_impulse_response_delay(impulse_response, N=1)[source]#

Find the delay in sub-sample values of an impulse response.

The method relies on the analytic part of the cross-correlation function of the impulse response and it’s minimum-phase equivalent, which is zero for the maximum of the correlation function. For sub-sample root finding, the analytic signal is approximated using a polynomial of order N. The algorithm is based on [3] with the following modifications:

Values with negative gradient used for polynomial fitting are rejected, allowing to use larger part of the signal for fitting.
By default a first order polynomial is used, as the slope of the analytic signal should in theory be linear.

Alternatively see pyfar.dsp.find_impulse_response_start. For complex-valued time signals, the delay is calculated separately for the real and complex part, and its minimum value returned.

Parameters:

impulse_response (Signal) – The impulse response.
N (int, optional) – The order of the polynom used for root finding, by default 1.

Returns:

delay – Delay of the impulse response, as an array of shape cshape. Can be floating point values in the case of sub-sample values.

Return type:

numpy.ndarray, float

References

Examples

Create a band-limited impulse shifted by 0.5 samples and estimate the starting sample of the impulse and plot.

>>> import pyfar as pf
>>> import numpy as np
>>> n_samples = 64
>>> delay_samples = n_samples // 2 + 1/2
>>> ir = pf.signals.impulse(n_samples)
>>> ir = pf.dsp.linear_phase(ir, delay_samples, unit='samples')
>>> start_samples = pf.dsp.find_impulse_response_delay(ir)
>>> ax = pf.plot.time(ir, unit='ms', label='impulse response')
>>> ax.axvline(
...     start_samples/ir.sampling_rate*1e3,
...     color='k', linestyle='-.', label='start sample')
>>> ax.legend()

(Source code, png, hires.png, pdf)

pyfar.dsp.find_impulse_response_start(impulse_response, threshold=20)[source]#

Find the start sample of an impulse response.

The start sample is identified as the first sample which is below the threshold level relative to the maximum level of the impulse response. For room impulse responses, ISO 3382 [4] specifies a threshold of 20 dB. This function is primary intended to be used when processing room impulse responses. For complex-valued time signals the onset is computed separately for the real and imaginary part. Alternatively see pyfar.dsp.find_impulse_response_delay.

Parameters:

impulse_response (pyfar.Signal) – The impulse response
threshold (float, optional) – The threshold level in dB, by default 20, which complies with ISO 3382.

Returns:

start_sample – Sample at which the impulse response starts

Return type:

numpy.ndarray, int

Notes

The function tries to estimate the PSNR in the IR based on the signal power in the last 10 percent of the IR. The automatic estimation may fail if the noise spectrum is not white or the impulse response contains non-linear distortions. If the PSNR is lower than the specified threshold, the function will issue a warning.

References

Examples

Create a band-limited impulse shifted by 0.5 samples and estimate the starting sample of the impulse and plot.

>>> import pyfar as pf
>>> import numpy as np
>>> n_samples = 256
>>> delay_samples = n_samples // 2 + 1/2
>>> ir = pf.signals.impulse(n_samples)
>>> ir = pf.dsp.linear_phase(ir, delay_samples, unit='samples')
>>> start_samples = pf.dsp.find_impulse_response_start(ir)
>>> ax = pf.plot.time(ir, unit='ms', label='impulse response', dB=True)
>>> ax.axvline(
...     start_samples/ir.sampling_rate*1e3,
...     color='k', linestyle='-.', label='start sample')
>>> ax.axhline(
...     20*np.log10(np.max(np.abs(ir.time)))-20,
...     color='k', linestyle=':', label='threshold')
>>> ax.legend()

(Source code, png, hires.png, pdf)

Create a train of weighted impulses with levels below and above the threshold, serving as a very abstract room impulse response. The starting sample is identified as the last sample below the threshold relative to the maximum of the impulse response.

>>> import pyfar as pf
>>> import numpy as np
>>> n_samples = 64
>>> delays = np.array([14, 22, 26, 30, 33])
>>> amplitudes = np.array([-35, -22, -6, 0, -9], dtype=float)
>>> ir = pf.signals.impulse(n_samples, delays, 10**(amplitudes/20))
>>> ir.time = np.sum(ir.time, axis=0)
>>> start_sample_est = pf.dsp.find_impulse_response_start(
...     ir, threshold=20)
>>> ax = pf.plot.time(
...     ir, dB=True, unit='samples',
...     label=f'peak samples: {delays}')
>>> ax.axvline(
...     start_sample_est, linestyle='-.', color='k',
...     label=f'ir start sample: {start_sample_est}')
>>> ax.axhline(
...     20*np.log10(np.max(np.abs(ir.time)))-20,
...     color='k', linestyle=':', label='threshold')
>>> ax.legend()

(Source code, png, hires.png, pdf)

pyfar.dsp.fractional_time_shift(signal, shift, unit='samples', order=30, side_lobe_suppression=60, mode='linear')[source]#

Apply fractional time shift to input data.

This function uses a windowed Sinc filter (Method FIR-2 in [5] according to Equations 21 and 22) to apply fractional delays, i.e., non-integer delays to an input signal. A Kaiser window according to [6] Equations (10.12) and (10.13) is used, which offers the possibility to control the side lobe suppression.

Parameters:

signal (Signal) – The input data
shift (float, array like) – The fractional shift in samples (positive or negative). If this is a float, the same shift is applied to all channels of signal. If this is an array like different delays are applied to the channels of signal. In this case it must broadcast to signal.cshape (see Numpy broadcasting)
unit (str, optional) – The unit of the shift. Either ‘samples’ or ‘s’. Defaults to ‘samples’.
order (int, optional) – The order of the fractional shift (sinc) filter. The precision of the filter increases with the order. High frequency errors decrease with increasing order. The order must be smaller than signal.n_samples. The default is 30.
side_lobe_suppression (float, optional) – The side lobe suppression of the Kaiser window in dB. The default is 60.
mode (str, optional) –
The filtering mode

"linear"
Apply linear shift, i.e., parts of the signal that are shifted to times smaller than 0 samples and larger than signal.n_samples disappear.

"cyclic"
Apply a cyclic shift, i.e., parts of the signal that are shifted to values smaller than 0 are wrapped around to the end, and parts that are shifted to values larger than signal.n_samples are wrapped around to the beginning.

The default is "linear"

Returns:

signal – The delayed input data

Return type:

Signal

References

Examples

Apply a fractional shift of 2.3 samples using filters of orders 6 and 30

>>> import pyfar as pf
>>> import matplotlib.pyplot as plt
>>>
>>> signal = pf.signals.impulse(64, 10)
>>>
>>> pf.plot.use()
>>> _, ax = plt.subplots(3, 1, figsize=(8, 8))
>>> pf.plot.time_freq(signal, ax=ax[:2], label="input", unit='ms')
>>> pf.plot.group_delay(signal, ax=ax[2], unit="samples")
>>>
>>> for order in [30, 6]:
>>>     delayed = pf.dsp.fractional_time_shift(
...         signal, 2.3, order=order)
>>>     pf.plot.time_freq(delayed, ax=ax[:2],
...                       label=f"delayed, order={order}", unit='ms')
>>>     pf.plot.group_delay(delayed, ax=ax[2], unit="samples")
>>>
>>> ax[1].set_ylim(-15, 5)
>>> ax[2].set_ylim(8, 14)
>>> ax[0].legend()

(Source code, png, hires.png, pdf)

Apply a shift that exceeds the signal length using the modes "linear" and "cyclic"

>>> import pyfar as pf
>>>
>>> signal = pf.signals.impulse(32, 16)
>>>
>>> ax = pf.plot.time(signal, label="input", unit='ms')
>>>
>>> for mode in ["cyclic", "linear"]:
>>>     delayed = pf.dsp.fractional_time_shift(
...         signal, 25.3, order=10, mode=mode)
>>>     pf.plot.time(delayed, label=f"delayed, mode={mode}", unit='ms')
>>>
>>> ax.legend()

(Source code, png, hires.png, pdf)

pyfar.dsp.group_delay(signal, frequencies=None, method='fft')[source]#

Returns the group delay of a signal in samples.

Parameters:

signal (Signal) – An audio signal object from the pyfar signal class
frequencies (array-like) – Frequency or frequencies in Hz at which the group delay is calculated. The default is None, in which case signal.frequencies is used.
method ('scipy', 'fft', optional) – Method to calculate the group delay of a Signal. Both methods calculate the group delay using the method presented in [7] avoiding issues due to discontinuities in the unwrapped phase. Note that the scipy version additionally allows to specify frequencies for which the group delay is evaluated. The default is 'fft', which is faster.

Returns:

group_delay – Frequency dependent group delay of shape (cshape, frequencies).

Return type:

numpy array

References

pyfar.dsp.kaiser_window_beta(A)[source]#

Return a shape parameter beta to create kaiser window based on desired side lobe suppression in dB.

This function can be used to call time_window with window=('kaiser', beta).

Parameters:: A (float) – Side lobe suppression in dB
Returns:: beta – Shape parameter beta after [8], Eq. 7.75
Return type:: float

References

pyfar.dsp.linear_phase(signal, group_delay, unit='samples')[source]#

Set the phase to a linear phase with a specified group delay.

The linear phase signal is computed as

\[H_{\mathrm{lin}} = |H| \mathrm{e}^{-j \omega \tau}\,,\]

with \(H\) the complex spectrum of the input data, \(|\cdot|\) the absolute values, \(\omega\) the frequency in radians and \(\tau\) the group delay in seconds.

Parameters:

signal (Signal) – input data
group_delay (float, array like) – The desired group delay of the linear phase signal according to unit. A reasonable value for most cases is signal.n_samples / 2 samples, which results in a time signal that is symmetric around the center. If group delay is a list or array it must broadcast with the channel layout of the signal (signal.cshape).
unit (string, optional) – Unit of the group delay. Can be 'samples' or 's' for seconds. The default is 'samples'.

Returns:

signal – linear phase copy of the input data

Return type:

Signal

pyfar.dsp.minimum_phase(signal, n_fft=None, truncate=True)[source]#

Calculate the minimum phase equivalent of a finite impulse response.

The method is based on the Hilbert transform of the real-valued cepstrum of the finite impulse response, that is the cepstrum of the magnitude spectrum only. As a result the magnitude spectrum is not distorted. Potential aliasing errors can occur due to the Fourier transform based calculation of the magnitude spectrum, which however are negligible if the length of Fourier transform n_fft is sufficiently high. [9] (Section 8.5.4)

Parameters:

signal (Signal) – The finite impulse response for which the minimum-phase version is computed.
n_fft (int, optional) – The FFT length used for calculating the cepstrum. Should be at least a few times larger than signal.n_samples. The default None uses eight times the signal length rounded up to the next power of two, that is: 2**int(np.ceil(np.log2(n_samples * 8))).
truncate (bool, optional) – If truncate is True, the resulting minimum phase impulse response is truncated to a length of signal.n_samples//2 + signal.n_samples % 2. This avoids aliasing described above in any case but might distort the magnitude response if signal.n_samples is to low. If truncate is False the output signal has the same length as the input signal. The default is True.

Returns:

signal_minphase – The minimum phase version of the input data.

Return type:

Signal

References

Examples

Create a minimum phase equivalent of a linear phase FIR low-pass filter

>>> import pyfar as pf
>>> import numpy as np
>>> from scipy.signal import remez
>>> import matplotlib.pyplot as plt
>>> freq = [0, 0.2, 0.3, 1.0]
>>> h_linear = pf.Signal(remez(151, freq, [1, 0], fs=2.), 44100)
>>> # create minimum phase impulse responses
>>> h_min = pf.dsp.minimum_phase(h_linear, truncate=False)
>>> # plot the result
>>> pf.plot.use()
>>> fig, axs = plt.subplots(3, figsize=(8, 6))
>>> pf.plot.time(h_linear, ax=axs[0], unit='ms')
>>> pf.plot.time(h_min, ax=axs[0], unit='ms')
>>> axs[0].grid(True)
>>> pf.plot.freq(h_linear, ax=axs[1])
>>> pf.plot.group_delay(h_linear, ax=axs[2], unit="ms")
>>> pf.plot.freq(h_min, ax=axs[1])
>>> pf.plot.group_delay(h_min, ax=axs[2], unit="ms")
>>> axs[2].legend(['Linear', 'Minimum'], loc=3, ncol=2)
>>> axs[2].set_ylim(-2.5, 2.5)

(Source code, png, hires.png, pdf)

pyfar.dsp.normalize(signal, reference_method='max', domain='auto', channel_handling='individual', target=1, limits=(None, None), unit=None, return_reference=False, nan_policy='raise')[source]#

Apply a normalization.

In the default case, the normalization ensures that the maximum absolute amplitude of the signal after normalization is 1. This is achieved by the multiplication

signal_normalized = signal * target / reference,

where target equals 1 and reference is the maximum absolute amplitude before the normalization.

Several normalizations are possible, which in fact are different ways of computing the reference value (e.g. based on the spectrum). See the parameters for details.

Parameters:

signal (Signal, TimeData, FrequencyData) – Input signal.
reference_method (string, optional) –
Reference method to compute the channel-wise reference value using the data according to domain.

'max'
Compute the maximum absolute value per channel.

'mean'
Compute the mean absolute values per channel.

'energy'
Compute the energy per channel using energy. Note that the square root of the energy is used as reference, after handling multi-channel signals (see channel_handling below). This is required for the energy of the normalized signal to match the target.

'power'
Compute the power per channel using power. Note that the square root of the power is used as reference, after handling multi-channel signals (see channel_handling below). This is required for the power of the normalized signal to match the target.

'rms'
Compute the RMS per channel using rms.

The default is 'max'.
domain (string) –
Determines which data is used to compute the reference value.

'time'
Use the absolute of the time domain data np.abs(signal.time).

'freq'
Use the magnitude spectrum np.abs(signal.freq). Note that the normalized magnitude spectrum is used pyfar examples gallery (cf. FFT normalization).

'auto'
Uses 'time' domain normalization for Signal and TimeData objects and 'freq' domain normalization for FrequencyData objects.

The default is 'auto'.
channel_handling (string, optional) –

Define how channel-wise reference values are handeled for multi-
channel signals. This parameter does not affect single-channel signals.

'individual'
Separate normalization of each channel individually.

'max'
Normalize to the maximum reference value across channels.

'min'
Normalize to the minimum reference value across channels.

'mean'
Normalize to the mean reference value across the channels.

The default is 'individual'.
target (scalar, array) – The target to which the signal is normalized. Can be a scalar or an array. In the latter case the shape of target must be broadcastable to cshape. The default is 1.
limits (tuple, array_like) – Restrict the time or frequency range that is used to compute the reference value. Two element tuple specifying upper and lower limit according to domain and unit. A None element means no upper or lower limitation. The default (None, None) uses the entire signal. Note that in case of limiting in samples or bins with unit=None, the second value defines the first sample/bin that is excluded. Also note that limits need to be (None, None) if reference_method is rms, power or energy.
unit (string, optional) –
Unit of limits.

's'
Set limits in seconds in case of time domain normalization. Uses find_nearest_time to find the limits.

'Hz'
Set limits in hertz in case of frequency domain normalization. Uses find_nearest_frequency

The default None assumes that limits is given in samples in case of time domain normalization and in bins in case of frequency domain normalization.
return_reference (bool) – If return_reference=True, the function also returns the reference values for the channels. The default is False.
nan_policy (string, optional) –
Define how to handle NaNs in input signal.

'propagate'
If the input signal includes NaNs within the time or frequency range , NaN will be used as normalization reference. The resulting output signal values are NaN.

'omit'
NaNs will be omitted in the normalization. Cshape will still remain, as the normalized signal still includes the NaNs.

'raise'
A ValueError will be raised, if the input signal includes NaNs.

The default is ‘raise’.

Returns:

normalized_signal (Signal, TimeData, FrequencyData) – The normalized input signal.
reference_norm (numpy.ndarray) – The reference values used for normalization. Only returned if return_reference is True.

Examples

Time domain normalization with default parameters

>>> import pyfar as pf
>>> signal = pf.signals.sine(1e3, 441, amplitude=2)
>>> signal_norm = pf.dsp.normalize(signal)
>>> # Plot input and normalized Signal
>>> ax = pf.plot.time(signal, label='Original Signal', unit='ms')
>>> pf.plot.time(signal_norm, label='Normalized Signal', unit='ms')
>>> ax.legend()

(Source code, png, hires.png, pdf)

Frequency normalization with a restricted frequency range and targed in dB

>>> import pyfar as pf
>>> sine1 = pf.signals.sine(1e3, 441, amplitude=2)
>>> sine2 = pf.signals.sine(5e2, 441, amplitude=.5)
>>> signal = sine1 + sine2
>>> # Normalize to dB target in restricted frequency range
>>> target_dB = 0
>>> signal_norm = pf.dsp.normalize(signal, target=10**(target_dB/20),
...     domain="freq", limits=(400, 600), unit="Hz")
>>> # Plot input and normalized Signal
>>> ax = pf.plot.time_freq(signal_norm, label='Normalized Signal',
...                        unit='ms')
>>> pf.plot.time_freq(signal, label='Original Signal', unit='ms')
>>> ax[1].set_ylim(-15, 15)
>>> ax[1].legend()

(Source code, png, hires.png, pdf)

pyfar.dsp.pad_zeros(signal, pad_width, mode='end')[source]#

Pad a signal with zeros in the time domain.

Parameters:

signal (Signal) – The signal which is to be extended.
pad_width (int) – The number of samples to be padded.
mode (str, optional) –
The padding mode:

'end'
Append zeros to the end of the signal

'beginning'
Prepend zeros to the beginning of the signal

'center'
Insert the number of zeros in the middle of the signal. This mode can be used to pad signals with a symmetry with respect to the time t=0.

The default is 'end'.

Returns:

The zero-padded signal.

Return type:

Signal

Examples

>>> import pyfar as pf
>>> impulse = pf.signals.impulse(512, amplitude=1)
>>> impulse_padded = pf.dsp.pad_zeros(impulse, 128, mode='end')

pyfar.dsp.phase(signal, deg=False, unwrap=False)[source]#

Returns the phase for a given signal object.

Parameters:

signal (Signal, FrequencyData) – pyfar Signal or FrequencyData object.
deg (Boolean) – Specifies, whether the phase is returned in degrees or radians.
unwrap (Boolean) – Specifies, whether the phase is unwrapped or not. If set to '360', the phase is wrapped to 2 pi.

Returns:

phase – The phase of the signal.

Return type:

numpy array

pyfar.dsp.power(signal)[source]#

Compute the power of a signal.

The power is calculated as

\[\frac{1}{N}\sum_{n=0}^{N-1}|x[n]|^2\]

based on the time data for each channel separately.

Parameters:: signal (Signal) – The signal to compute the power from.
Returns:: data – The channel-wise power of the input signal.
Return type:: numpy.ndarray

Notes

Due to the calculation based on the time data, the returned power is independent of the signal’s fft_norm. The power equals the squared RMS of a signal. energy and rms can be used to compute the energy and the RMS.

pyfar.dsp.regularized_spectrum_inversion(signal, frequency_range, regu_outside=1.0, regu_inside=1e-10, regu_final=None, normalized=True)[source]#

Invert the spectrum of a signal applying frequency dependent regularization.

Regularization can either be specified within a given frequency range using two different regularization factors, or for each frequency individually using the parameter regu_final. In the first case the regularization factors for the frequency regions are cross-faded using a raised cosine window function with a width of \(\sqrt{2}f\) above and below the given frequency range. Note that the resulting regularization function is adjusted to the quadratic maximum of the given signal. In case the regu_final parameter is used, all remaining options are ignored and an array matching the number of frequency bins of the signal needs to be given. In this case, no normalization of the regularization function is applied.

Finally, the inverse spectrum is calculated as [10], [11],

\[S^{-1}(f) = \frac{S^*(f)}{S^*(f)S(f) + \epsilon(f)}\]

Parameters:

signal (Signal) – The signals which spectra are to be inverted.
frequency_range (tuple, array_like, double) – The upper and lower frequency limits outside of which the regularization factor is to be applied.
regu_outside (float, optional) – The normalized regularization factor outside the frequency range. The default is 1.
regu_inside (float, optional) – The normalized regularization factor inside the frequency range. The default is 10**(-200/20) (-200 dB).
regu_final (float, array_like, optional) – The final regularization factor for each frequency, default None. If this parameter is set, the remaining regularization factors are ignored.
normalized (bool) – Flag to indicate if the normalized spectrum (according to signal.fft_norm) should be inverted. The default is True.
freq_range (tuple, array_like, double) – The upper and lower frequency limits outside of which the regularization factor is to be applied. 'freq_range' parameter will be deprecated in pyfar 0.8.0 in favor of 'frequency_range'.

Returns:

The resulting signal after inversion.

Return type:

Signal

References

pyfar.dsp.resample(signal, sampling_rate, match_amplitude='auto', frac_limit=None, post_filter=False)[source]#

Resample signal to new sampling rate.

The SciPy function scipy.signal.resample_poly is used for resampling. The resampling ratio L = sampling_rate/signal.sampling_rate is approximated by a fraction of two integer numbers up/down to first upsample the signal by up and then downsample by down. This way up and down are smaller than the respective new and old sampling rates.

Note

sampling_rate should be divisible by 10, otherwise it can cause an infinite loop in the resample_poly function.

The amplitudes of the resampled signal can match the amplitude of the input signal in the time or frequency domain. See the parameter match_amplitude and the examples for more information.

Parameters:

signal (Signal) – Input data to be resampled
sampling_rate (number) – The new sampling rate in Hz
match_amplitude (string) –
Define the domain to match the amplitude of the resampled data.

'auto'
Chooses domain to maintain the amplitude automatically, depending on the signal.signal_type. Sets match_amplitude == 'freq' for signal.signal_type = 'energy' like impulse responses and match_amplitude == 'time' for other signals.

'time'
Maintains the amplitude in the time domain. This is useful for recordings such as speech or music and must be used if signal.signal_type = 'power'.

'freq'
Maintains the amplitude in the frequency domain by multiplying the resampled signal by 1/L (see above). This is often desired when resampling impulse responses.

The default is 'auto'.
frac_limit (int) –
Limit the denominator for approximating the resampling factor L (see above). This can be used in case the resampling gets stuck in an infinite loop (see note above) at the potenital cost of not exactly realizing the target sampling rate.

The default is None, which uses frac_limit = 1e6.
post_filter (bool, optional) – In some cases the up-sampling causes artifacts above the Nyquist frequency of the input signal, i.e., signal.sampling_rate/2. If True the artifacts are suppressed by applying a zero-phase Elliptic filter with a pass band ripple of 0.1 dB, a stop band attenuation of 60 dB. The pass band edge frequency is signal.sampling_rate/2. The stop band edge frequency is the minimum of 1.05 times the pass band frequency and the new Nyquist frequency (sampling_rate/2). The default is False. Note that this is only applied in case of up-sampling.

Returns:

signal – The resampled signal of the input data with a length of up/down * signal.n_samples samples.

Return type:

pyfar.Signal

Examples

For power signals, the amplitude of the resampled signal is automatically correct in the time and frequency domain if match_amplitude="time"

>>> import pyfar as pf
>>> import matplotlib.pyplot as plt
>>>
>>> signal = pf.signals.sine(200, 4800, sampling_rate=48000)
>>> resampled = pf.dsp.resample(signal, 96000)
>>>
>>> pf.plot.time_freq(signal, label="original", unit='ms')
>>> pf.plot.time_freq(resampled, c="y", ls=":", unit='ms',
...                   label="resampled (time domain matched)")
>>> plt.legend()

(Source code, png, hires.png, pdf)

With some energy signals, such as impulse responses, the amplitude can only be correct in the time or frequency domain due to the lack of normalization by the number of samples. In such cases, it is often desired to match the amplitude in the frequency domain

>>> import pyfar as pf
>>> import matplotlib.pyplot as plt
>>>
>>> signal = pf.signals.impulse(128, 64, sampling_rate=48000)
>>> resampled_time = pf.dsp.resample(
...     signal, 96000, match_amplitude = "time")
>>> resampled_freq = pf.dsp.resample(
...     signal, 96000, match_amplitude = "freq")
>>>
>>> pf.plot.time_freq(signal, label="original", unit='ms')
>>> pf.plot.time_freq(resampled_freq, dashes=[2, 3], unit='ms',
...                   label="resampled (freq. domain matched)")
>>> ax = pf.plot.time_freq(resampled_time, ls=":", unit='ms',
...                   label="resampled (time domain matched)", c='y')
>>> ax[0].set_xlim(1.2,1.46)
>>> plt.legend()

(Source code, png, hires.png, pdf)

pyfar.dsp.rms(signal)[source]#

Compute the root mean square (RMS) of a signal.

The RMS is calculated as

\[\sqrt{\frac{1}{N}\sum_{n=0}^{N-1}|x[n]|^2}\]

based on the time data for each channel separately.

Parameters:: signal (Signal) – The signal to compute the RMS from.
Returns:: data – The channel-wise RMS of the input signal.
Return type:: numpy.ndarray

Notes

The RMS equals the square root of the signal’s power. energy and power can be used to compute the energy and the power.

pyfar.dsp.smooth_fractional_octave(signal, num_fractions, mode='magnitude_zerophase', window='boxcar')[source]#

Smooth spectrum with a fractional octave width.

The smoothing is done according to Tylka et al. 2017 [12] (method 2) in three steps:

Interpolate the spectrum to a logarithmically spaced frequency scale
Smooth the spectrum by convolution with a smoothing window
Interpolate the spectrum to the original linear frequency scale

Smoothing of complex-valued time data is not implemented.

Parameters:

signal (pyfar.Signal) – The input data.
num_fractions (number) – The width of the smoothing window in fractional octaves, e.g., 3 will apply third octave smoothing and 1 will apply octave smoothing.
mode (str, optional) –

"magnitude_zerophase"
Only the magnitude response, i.e., the absolute spectrum is smoothed. Note that this return a zero-phase signal. It might be necessary to generate a minimum or linear phase if the data is subject to further processing after the smoothing (cf. minimum_phase and linear_phase)

"magnitude"
Smooth the magnitude and keep the phase of the input signal.

"magnitude_phase"
Separately smooth the magnitude and unwrapped phase response.

"complex"
Separately smooth the real and imaginary part of the spectrum.

Note that the modes magnitude_zerophase and magnitude make sure that the smoothed magnitude response is as expected at the cost of an artificial phase response. This is often desired, e.g., when plotting signals or designing compensation filters. The modes magnitude_phase and complex smooth all information but might cause a high frequency energy loss in the smoothed magnitude response. The default is "magnitude_zerophase".
window (str, optional) – String that defines the smoothing window. All windows from time_window that do not require an additional parameter can be used. The default is “boxcar”, which uses the most commonly used rectangular window.

Returns:

signal (pyfar.Signal) – The smoothed output data
window_stats (tuple) – A tuple containing information about the smoothing process

n_window
The window length in (logarithmically spaced) samples

num_fractions
The actual width of the window in fractional octaves. This can deviate from the desired width because the smoothing window must have an integer sample length

Notes

Method 3 in Tylka at al. 2017 is mathematically more elegant at the price of a largely increased computational and memory cost. In most practical cases, methods 2 and 3 yield close to identical results (cf. Fig. 2 and 3 in Tylka et al. 2017). If the spectrum contains extreme discontinuities, however, method 3 is superior (see examples below).

References

Examples

Octave smoothing of continuous spectrum consisting of two bell filters.

>>> import pyfar as pf
>>> signal = pf.signals.impulse(441)
>>> signal = pf.dsp.filter.bell(signal, 1e3, 12, 1, "III")
>>> signal = pf.dsp.filter.bell(signal, 10e3, -60, 100, "III")
>>> smoothed, _ = pf.dsp.smooth_fractional_octave(signal, 1)
>>> ax = pf.plot.freq(signal, label="input")
>>> pf.plot.freq(smoothed, label="smoothed")
>>> ax.legend(loc=3)

(Source code, png, hires.png, pdf)

Octave smoothing of the discontinuous spectrum of a sine signal causes artifacts at the edges due to the intermediate interpolation steps (cf. Tylka et al. 2017, Fig. 4). However this is a rather unusual application and is mentioned only for the sake of completeness.

>>> import pyfar as pf
>>> signal = pf.signals.sine(1e3, 4410)
>>> signal.fft_norm = "amplitude"
>>> smoothed, _ = pf.dsp.smooth_fractional_octave(signal, 1)
>>> ax = pf.plot.freq(signal, label="input")
>>> pf.plot.freq(smoothed, label="smoothed")
>>> ax.set_xlim(200, 4e3)
>>> ax.set_ylim(-45, 5)
>>> ax.legend(loc=3)

(Source code, png, hires.png, pdf)

pyfar.dsp.soft_limit_spectrum(signal, limit, knee, frequency_range=None, direction='upper', log_prefix=None)[source]#

Soft limiting the magniude spectrum.

Soft limiting gradually increases the gain reduction to avoid discontinuities in the data that would appear in hard limiting. The transition between the magnitude where no limiting is applied to the magnitude where the limiting reaches its full effect is termed knee (see examples below).

Note that the limiting is applied on signal.freq, i.e. the data after the FFT normalization.

Parameters:

signal (Signal, FrequencyData) – The input data
limit (number, array like) – The gain in dB at which the limiting reaches its full effect. If this is a number, the same limit is applied to all frequencies. If this an array like, it must be broadcastable to signal.freq to apply frequency-dependent limits.
knee (number, string) – If this is a number, a knee with a width of number dB according to [13] Eq. (4) is applied. This definition of the knee originates from the classic limiting audio effect. If this is 'arctan' an arcus tangens knee according to [14] Section 3.6.4 is applied. This knee definition originates from microphone array signal processing.
frequency_range (array like, optional) – Frequency range in which the limiting is applied. This must be an array like containing the lower and upper limit in Hz. The default None applies the limiting to all frequencies.
direction (str, optional) –
Define how the limiting works

'upper' (default)
Soft limiting signal to enforce an aboslute maximum value of limit.

'lower'
Soft limiting signal to enforce an absolute minimum value of limit
log_prefix (float, int) – The log prefix is used to linearize the limit and knee, e.g., limit_linear = 10**(limit / log_prefix).The default None, uses 10 for signals with 'psd' and 'power' FFT normalization and 20 otherwise.

Returns:

limited – The limited copy of the input data.

Return type:

Signal, FrequencyData

Examples

Illustrate effect of limit and knee

>>> import pyfar as pf
>>> import numpy as np
>>>
>>> signal = pf.FrequencyData(
...     10**(np.arange(-20, 21)/20), np.arange(0, 41))
>>> for knee in ['arctan', 20, 0]:
>>>     limited = pf.dsp.soft_limit_spectrum(signal, 0, knee)
>>>     pf.plot.freq(limited, freq_scale='linear', label=f'{knee=}')
>>> ax = pf.plot.freq(signal, freq_scale='linear',
...                   linestyle='--', label='input')
>>>
>>> ax.legend(loc='upper left')
>>> ax.set_xlim(0, 40)
>>> ax.set_ylim(-20, 20)
>>> bbox=dict(facecolor='white', edgecolor='none', alpha=0.5)
>>> arrowprops= dict(
...     arrowstyle="<->", shrinkA=0, shrinkB=0,
...     color=pf.plot.color('r'))
>>> ax.annotate("limit=0 dB", (20, 0), (10, 0), bbox=bbox,
...             arrowprops=dict(arrowstyle="->"), va='center')
>>> ax.text(20, -14, 'knee=20 dB', ha='center',
...         color=pf.plot.color('r'), bbox=bbox)
>>> ax.annotate("", (10, -15), (30, -15), arrowprops=arrowprops)

(Source code, png, hires.png, pdf)

Apply limiting with knee (soft limiting) and without knee (hard limiting)

>>> import pyfar as pf
>>>
>>> signal = pf.signals.impulse(4096)
>>> signal = pf.dsp.filter.bell(signal, 1e3, 20, 1)
>>>
>>> soft = pf.dsp.soft_limit_spectrum(signal, limit=10, knee=12)
>>> hard = pf.dsp.soft_limit_spectrum(signal, limit=10, knee=0)
>>>
>>> pf.plot.freq(soft, label='soft (knee > 0)')
>>> pf.plot.freq(hard, label='hard (knee = 0)')
>>> ax = pf.plot.freq(signal, ls='--', label='original')
>>>
>>> ax.set_ylim(-5, 25)
>>> ax.legend()

(Source code, png, hires.png, pdf)

Apply soft limiting above and below +/-10 dB

>>> import pyfar as pf
>>>
>>> signal = pf.signals.impulse(4096)
>>> signal = pf.dsp.filter.bell(signal, 2e2, -20, 3)
>>> signal = pf.dsp.filter.bell(signal, 6e3, 20, 3)
>>>
>>> upper = pf.dsp.soft_limit_spectrum(
...     signal, limit=10, knee=6, direction="upper")
>>> lower = pf.dsp.soft_limit_spectrum(
...     signal, limit=-10, knee=6, direction="lower")
>>>
>>> pf.plot.freq(upper, label='upper limit')
>>> pf.plot.freq(lower, label='lower limit')
>>> ax = pf.plot.freq(signal, ls='--', label='original')
>>>
>>> ax.set_ylim(-25, 25)
>>> ax.legend()

(Source code, png, hires.png, pdf)

Apply frequency dependent soft limiting

>>> import pyfar as pf
>>>
>>> # input data
>>> signal = pf.signals.impulse(2**13)
>>> signal = pf.dsp.filter.bell(signal, 2e2, 20, 4)
>>> signal = pf.dsp.filter.bell(signal, 2e3, 20, 4)
>>>
>>> # frequency dependent limit
>>> limit = pf.dsp.filter.low_shelf(
...     pf.signals.impulse(2**13), 500, -10, 2) * 10
>>> limit = pf.dsp.decibel(limit).flatten()
>>>
>>> # soft limiting input signal
>>> limited = pf.dsp.soft_limit_spectrum(signal, limit, 10)
>>>
>>> # plot
>>> ax = pf.plot.freq(signal, label='original')
>>> pf.plot.freq(limited, label='limited')
>>> ax.plot(signal.frequencies, limit, label='limit', c='k', ls='--')
>>> ax.set_ylim(-5, 25)
>>> ax.legend(loc='upper left')

(Source code, png, hires.png, pdf)

References

pyfar.dsp.spectrogram(signal, window='hann', window_length=1024, window_overlap_fct=0.5, normalize=True)[source]#

Compute the magnitude spectrum versus time.

This is a wrapper for scipy.signal.ShortTimeFFT.spectrogram with two differences. First, the returned times refer to the start of the FFT blocks, i.e., the first time is always 0 whereas it is window_length/2 in scipy. Second, the returned spectrogram is normalized according to signal.fft_norm if the normalize parameter is set to True.

Parameters:

signal (Signal) – Signal to compute spectrogram of.
window (str) – Specifies the window (see scipy.signal.windows). The default is 'hann'.
window_length (integer) – Window length in samples, the default ist 1024.
window_overlap_fct (double) – Ratio of points to overlap between FFT segments [0…1]. The default is 0.5.
normalize (bool) – Flag to indicate if the FFT normalization should be applied to the spectrogram according to signal.fft_norm. The default is True.

Returns:

frequencies (numpy array) – Frequencies in Hz at which the magnitude spectrum was computed. For complex valued signals, frequencies and spectrogram is arranged such that 0 Hz bin is centered.
times (numpy array) – Times in seconds at which the magnitude spectrum was computed
spectrogram (numpy array)

pyfar.dsp.time_shift(signal, shift, mode='cyclic', unit='samples', pad_value=0.0)[source]#

Apply a cyclic or linear time-shift to a signal.

This function only allows integer value sample shifts. If unit 'time' is used, the shift samples will be rounded to the nearest integer value. For a shift using fractional sample values see fractional_time_shift.

Parameters:

signal (Signal) – The signal to be shifted
shift (int, float, array_like) – The time-shift value. A positive value will result in right shift on the time axis (delaying of the signal), whereas a negative value yields a left shift on the time axis (non-causal shift to a earlier time). If a single value is given, the same time shift will be applied to each channel of the signal. Individual time shifts for each channel can be performed by passing an array broadcastable to the signals channel dimensions cshape.
mode (str, optional) –
The shifting mode

"linear"
Apply linear shift, i.e., parts of the signal that are shifted to times smaller than 0 samples and larger than signal.n_samples disappear. To maintain the shape of the signal, the signal is padded at the respective other end. The pad value is determined by pad_type.

"cyclic"
Apply a cyclic shift, i.e., parts of the signal that are shifted to values smaller than 0 are wrapped around to the end, and parts that are shifted to values larger than signal.n_samples are wrapped around to the beginning.

The default is "cyclic"
unit (str, optional) – Unit of the shift variable, this can be either 'samples' or 's' for seconds. By default 'samples' is used. Note that in the case of specifying the shift time in seconds, the value is rounded to the next integer sample value to perform the shift.
pad_value (numeric, optional) – The pad value for linear shifts, by default 0. is used. Pad numpy.nan to the respective channels if the rms value of the signal is to be maintained for block-wise rms estimation of the noise power of a signal. Note that if NaNs are padded, the returned data will be a TimeData instead of Signal object.

Returns:

The time-shifted signal. This is a TimeData object in case a linear shift was done and the signal was padded with Nans. In all other cases, a Signal object is returned.

Return type:

Signal, TimeData

Examples

Individually do a cyclic shift of a set of ideal impulses stored in three different channels and plot the resulting signals

>>> import pyfar as pf
>>> import matplotlib.pyplot as plt
>>> # generate and shift the impulses
>>> impulse = pf.signals.impulse(
...     32, amplitude=(1, 1.5, 1), delay=(14, 15, 16))
>>> shifted = pf.dsp.time_shift(impulse, [-2, 0, 2])
>>> # time domain plot
>>> pf.plot.use('light')
>>> _, axs = plt.subplots(2, 1)
>>> pf.plot.time(impulse, ax=axs[0], unit='samples')
>>> pf.plot.time(shifted, ax=axs[1], unit='samples')
>>> axs[0].set_title('Original signals')
>>> axs[1].set_title('Shifted signals')

(Source code, png, hires.png, pdf)

Perform a linear time shift instead and pad with NaNs

>>> import pyfar as pf
>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> # generate and shift the impulses
>>> impulse = pf.signals.impulse(
...     32, amplitude=(1, 1.5, 1), delay=(14, 15, 16))
>>> shifted = pf.dsp.time_shift(
...     impulse, [-2, 0, 2], mode='linear', pad_value=np.nan)
>>> # time domain plot
>>> pf.plot.use('light')
>>> _, axs = plt.subplots(2, 1)
>>> pf.plot.time(impulse, ax=axs[0], unit='samples')
>>> pf.plot.time(shifted, ax=axs[1], unit='samples')
>>> axs[0].set_title('Original signals')
>>> axs[1].set_title('Shifted signals')

(Source code, png, hires.png, pdf)

pyfar.dsp.time_window(signal, interval, window='hann', shape='symmetric', unit='samples', crop='none', return_window=False)[source]#

Apply time window to signal.

This function uses the windows implemented in scipy.signal.windows.

Parameters:

signal (Signal) – Signal object to be windowed.
interval (array_like) – If interval has two entries, these specify the beginning and the end of the symmetric window or the fade-in / fade-out (see parameter shape). If interval has four entries, a window with fade-in between the first two entries and a fade-out between the last two is created, while it is constant in between (ignores shape). The unit of interval is specified by the parameter unit. See below for more details.
window (string, float, or tuple, optional) – The type of the window. See below for a list of implemented windows. The default is 'hann'.
shape (string, optional) –

'symmetric'
General symmetric window, the two values in interval define the first and last samples of the window.

'symmetric_zero'
Symmetric window with respect to t=0, the two values in interval define the first and last samples of fade-out. crop is ignored.

'left'
Fade-in, the beginning and the end of the fade is defined by the two values in interval. See Notes for more details.

'right'
Fade-out, the beginning and the end of the fade is defined by the two values in interval. See Notes for more details.

The default is 'symmetric'.
unit (string, optional) – Unit of interval. Can be set to 'samples' or 's' (seconds). Time values are rounded to the nearest sample. The default is 'samples'.
crop (string, optional) –

'none'
The length of the windowed signal stays the same.

'window'
The signal is truncated to the windowed part.

'end'
Only the zeros at the end of the windowed signal are cropped, so the original phase is preserved.

The default is 'none'.
return_window (bool, optional) – If True, both the windowed signal and the time window are returned. The default is False.

Returns:

signal_windowed (Signal) – Windowed signal object
window (Signal) – Time window used to create the windowed signal, only returned if return_window=True.

Notes

For a fade-in, the indexes of the samples given in interval denote the first sample of the window which is non-zero and the first which is one. For a fade-out, the samples given in interval denote the last sample which is one and the last which is non-zero.

This function calls scipy.signal.windows.get_window to create the window. Available window types:

boxcar
triang
blackman
hamming
hann
bartlett
flattop
parzen
bohman
blackmanharris
nuttall
barthann
kaiser (needs beta, see kaiser_window_beta)
gaussian (needs standard deviation)
general_gaussian (needs power, width)
dpss (needs normalized half-bandwidth)
chebwin (needs attenuation)
exponential (needs center, decay scale)
tukey (needs taper fraction)
taylor (needs number of constant sidelobes, sidelobe level)

If the window requires no parameters, then window can be a string. If the window requires parameters, then window must be a tuple with the first argument the string name of the window, and the next arguments the needed parameters.

Examples

Options for parameter shape.

>>> import pyfar as pf
>>> import numpy as np
>>> signal = pf.Signal(np.ones(100), 44100)
>>> for shape in ['symmetric', 'symmetric_zero', 'left', 'right']:
>>>     signal_windowed = pf.dsp.time_window(
...         signal, interval=[25,45], shape=shape)
>>>     ax = pf.plot.time(signal_windowed, label=shape, unit='ms')
>>> ax.legend(loc='right')

(Source code, png, hires.png, pdf)

Window with fade-in and fade-out defined by four values in interval.

>>> import pyfar as pf
>>> import numpy as np
>>> signal = pf.Signal(np.ones(100), 44100)
>>> signal_windowed = pf.dsp.time_window(
...         signal, interval=[25, 40, 60, 90], window='hann')
>>> pf.plot.time(signal_windowed, unit='ms')

(Source code, png, hires.png, pdf)

pyfar.dsp.wrap_to_2pi(x)[source]#

Wraps phase to 2 pi.

Parameters:: x (double) – Input phase to be wrapped to 2 pi.
Returns:: x – Phase wrapped to 2 pi`.
Return type:: double

pyfar.dsp.zero_phase(signal)[source]#

Calculate zero phase signal.

The zero phase signal is obtained by taking the absolute values of the spectrum

\[H_z = |H| = \sqrt{\mathrm{real}(H)^2 + \mathrm{imag}(H)^2},\]

where \(H\) is the complex valued spectrum of the input data and \(H_z\) the real valued zero phase spectrum.

The time domain data of a zero phase signal is symmetric around the first sample, e.g., signal.time[0, 1] == signal.time[0, -1].

Parameters:: signal (Signal, FrequencyData) – input data
Returns:: signal – zero phase copy of the input data
Return type:: Signal, FrequencyData