engineeringkit.signal_processing#

Domain hub · History · Examples · API reference

Sampling, aliasing and sinc reconstruction; uniform quantization and SNR; decimation and interpolation; the DFT, a hand-written radix-2 FFT, Fourier partial sums and the DCT; linear and circular convolution, overlap-add and correlation; the periodogram, windowed spectra, Welch’s averaged periodogram, the short-time Fourier transform, Goertzel’s algorithm and the cepstrum; autoregressive models by the Levinson-Durbin recursion; the Hilbert transform and analytic signal; Butterworth, Chebyshev and elliptic filters as polynomials or second-order sections (scipy.signal) with minimum-order formulas, the bilinear transform with prewarping, window-method FIR filters with Kaiser’s formulas and Parks-McClellan equiripple filters (scipy.signal.remez); the FIR Wiener filter (scipy.linalg.solve_toeplitz) and the LMS adaptive filter; and the Haar wavelet transform.

engineeringkit.signal_processing: sampling, Fourier analysis, spectra, filters, adaptive filters and wavelets.

Sampling, aliasing and sinc reconstruction; uniform quantization and SNR; decimation and interpolation; the DFT, a hand-written radix-2 FFT, Fourier partial sums and the DCT; linear and circular convolution, overlap-add and correlation; the periodogram, windowed spectra, Welch’s averaged periodogram, the short-time Fourier transform, Goertzel’s algorithm and the cepstrum; autoregressive models by the Levinson-Durbin recursion; the Hilbert transform and analytic signal; Butterworth, Chebyshev and elliptic IIR filters as polynomials or second-order sections, the bilinear transform with prewarping, window-method FIR filters with Kaiser’s design formulas and Parks-McClellan equiripple filters; the FIR Wiener filter and the LMS adaptive filter; and the Haar wavelet transform. NumPy’s FFT and SciPy’s filter design and filtering are used where they are the standard tool; the algorithms named after their inventors are written out.

class engineeringkit.signal_processing.ARModel(coefficients, noise_variance, reflection_coefficients, prediction_errors)[source]#

Bases: object

Autoregressive model \(x_n + \sum_{k=1}^{p} a_k x_{n-k} = e_n\).

Returned by engineeringkit.signal_processing.levinson_durbin() and engineeringkit.signal_processing.yule_walker().

Parameters:

  • coefficients (ndarray)

  • noise_variance (float)

  • reflection_coefficients (ndarray)

  • prediction_errors (ndarray)

coefficients: ndarray#

\([1, a_1, \ldots, a_p]\), the prediction-error filter \(A(z)\) [-].

Type:

ndarray

noise_variance: float#

Variance \(\sigma^2 = E_p\) of the white input \(e_n\) [signal units squared].

Type:

float

property order: int#

Model order \(p\).

Type:

int

plot(ax=None, fs=1.0)[source]#

Model spectrum in decibels against frequency; returns (fig, ax).

prediction_errors: ndarray#

Mean squared prediction error \(E_0..E_p\) of each order [signal units squared].

Type:

ndarray

psd(fs=1.0, n=512)[source]#

One-sided power spectral density \(P(f) = 2\sigma^2/(f_s|A(e^{j2\pi f/f_s})|^2)\).

Parameters:

  • fs (float, optional) – Sampling rate [Hz], positive.

  • n (int, optional) – Number of frequencies from 0 to \(f_s/2\), at least 2.

Returns:

(frequency, psd) [Hz] and [signal units squared per Hz]; its integral over \([0, f_s/2]\) is the model’s variance.

Return type:

tuple of ndarray

reflection_coefficients: ndarray#

\(k_1..k_p\), all inside (-1, 1) for a stable model [-].

Type:

ndarray

summary()[source]#

Order, noise variance, prediction gain \(E_0/E_p\) and whether every \(|k_m| < 1\).

Return type:

dict

class engineeringkit.signal_processing.LMSResult(output, error, weights, weight_history)[source]#

Bases: object

Output of engineeringkit.signal_processing.lms_filter().

Parameters:

  • output (ndarray)

  • error (ndarray)

  • weights (ndarray)

  • weight_history (ndarray)

error: ndarray#

Error \(e_n = d_n - y_n\) [signal units].

Type:

ndarray

output: ndarray#

Filter output \(y_n = w_n^T u_n\) [signal units].

Type:

ndarray

plot(ax=None)[source]#

Learning curve, the squared error against sample index on a log axis; returns (fig, ax).

summary()[source]#

Number of taps and mean squared error over the first and last tenth of the run.

Return type:

dict

weight_history: ndarray#

Weights before each update, shape (n_samples, n_taps) [-].

Type:

ndarray

weights: ndarray#

Weight vector after the last sample [-].

Type:

ndarray

class engineeringkit.signal_processing.STFTResult(time, frequency, values)[source]#

Bases: object

Output of engineeringkit.signal_processing.stft().

Parameters:

  • time (ndarray)

  • frequency (ndarray)

  • values (ndarray)

frequency: ndarray#

One-sided bin frequencies [Hz].

Type:

ndarray

plot(ax=None)[source]#

Spectrogram \(|X|^2\) in decibels relative to its peak; returns (fig, ax).

summary()[source]#

Frame and bin counts, and the time step and frequency resolution of the grid.

Return type:

dict

time: ndarray#

Centre time of each frame [s].

Type:

ndarray

values: ndarray#

Coefficients, shape (frequency.size, time.size) [signal units].

Type:

ndarray of complex

class engineeringkit.signal_processing.Spectrum(frequency, magnitude)[source]#

Bases: object

One-sided, unnormalized FFT magnitude of a real signal.

Parameters:

  • frequency (ndarray)

  • magnitude (ndarray)

frequency: ndarray#

Bin frequencies [Hz], from 0 to the Nyquist frequency.

Type:

ndarray

classmethod from_signal(x, fs)[source]#

Magnitude spectrum of real samples.

Parameters:

  • x (array_like) – Real, finite samples.

  • fs (float) – Sampling rate [Hz], positive.

Returns:

One-sided frequencies [Hz] and unnormalized FFT magnitudes.

Return type:

Spectrum

magnitude: ndarray#

\(|X_k|\) of the unnormalized real FFT [signal units].

Type:

ndarray

plot(ax=None)[source]#

Magnitude against frequency; returns (fig, ax).

summary()[source]#

Number of bins, bin spacing [Hz] and frequency of the largest bin [Hz].

Return type:

dict

engineeringkit.signal_processing.alias_frequency(f, fs)[source]#

Apparent frequency in \([0, f_s/2]\) of a tone at f sampled at fs.

Parameters:

  • f (float or array_like) – Tone frequency [Hz]; negative frequencies fold like positive ones.

  • fs (float) – Sampling rate [Hz], positive.

Returns:

Alias frequency \(|((f + f_s/2) \bmod f_s) - f_s/2|\) [Hz].

Return type:

float or ndarray

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> float(alias_frequency(70, 100))
30.0
engineeringkit.signal_processing.analytic_signal(x)[source]#

Analytic signal \(z = x + j\mathcal{H}x\) by suppressing the negative-frequency half of the DFT.

Parameters:

x (array_like) – Real 1D samples [signal units], at least 2.

Returns:

\(z\), with \(\operatorname{Re} z = x\) [signal units].

Return type:

ndarray of complex

Notes

The DFT bins are multiplied by 1 at DC (and at Nyquist for even \(N\)), 2 at positive frequencies and 0 at negative ones, then inverted. The record is treated as one period of a periodic signal, so a record that does not hold a whole number of cycles leaks at its ends. Equal to scipy.signal.hilbert. See Signal processing: sampled energy, spectra and quantization.

References

D. Gabor, “Theory of communication,” J. IEE 93, part III (1946) 429-457; S. L. Marple Jr., “Computing the discrete-time ‘analytic’ signal via FFT,” IEEE Trans. Signal Process. 47 (1999) 2600-2603.

Examples

>>> z = analytic_signal(np.cos(2 * np.pi * np.arange(8) / 8))
>>> (np.round(z.imag, 12) + 0.0).tolist()[:3]
[0.0, 0.707106781187, 1.0]
engineeringkit.signal_processing.bilinear_transform(b, a, fs, prewarp=None)[source]#

Digital filter \(H(z)\) from an analog \(H(s)\) by the bilinear transform, written out.

Parameters:

  • b (array_like) – Analog numerator and denominator coefficients in descending powers of \(s\) [-], with a of degree at least that of b.

  • a (array_like) – Analog numerator and denominator coefficients in descending powers of \(s\) [-], with a of degree at least that of b.

  • fs (float) – Sampling rate [Hz], positive.

  • prewarp (float, optional) – Frequency [Hz], below \(f_s/2\), at which the digital response is made to equal the analog one exactly. Without it \(K = 2f_s\).

Returns:

(b, a): coefficients in \(z^{-1}\) [-], with a[0] = 1.

Return type:

tuple of ndarray

Notes

Substituting \(s = K(1 - z^{-1})/(1 + z^{-1})\) and multiplying through by \((1 + z^{-1})^N\), each term \(c_m s^m\) of degree \(N\) polynomials becomes \(c_m K^m (1 - z^{-1})^m (1 + z^{-1})^{N-m}\). The left half \(s\)-plane maps inside the unit circle, so stable analog filters give stable digital ones, and the \(j\Omega\) axis maps onto the circle with \(\Omega = K\tan(\omega/2)\). Prewarping at \(f_0\) sets \(K = 2\pi f_0/\tan(\pi f_0/f_s)\). The result equals scipy.signal.bilinear with fs = K/2. See prewarp_frequency() and Signal processing: sampled energy, spectra and quantization.

References

R. M. Golden and J. F. Kaiser, “Design of wideband sampled-data filters,” Bell Syst. Tech. J. 43 (1964) 1533-1546.

Examples

>>> b, a = bilinear_transform([1.0], [1.0, 1.0], 0.5)
>>> b.tolist(), a.tolist()
([0.5, 0.5], [1.0, 0.0])
engineeringkit.signal_processing.butterworth_highpass(order, cutoff, fs, output='ba')[source]#

Digital Butterworth high-pass filter by the bilinear transform.

Parameters:

  • order (int) – Filter order \(n\), at least 1.

  • cutoff (float) – -3 dB frequency \(f_c\) [Hz], between 0 and \(f_s/2\).

  • fs (float) – Sampling rate [Hz], positive.

  • output ({"ba", "sos"}, optional) – Polynomials or second-order sections, as in butterworth_lowpass().

Returns:

(b, a) or second-order sections [-].

Return type:

tuple of ndarray or ndarray

Notes

\(|H(f)|^2 = 1/(1 + (\tan(\pi f_c/f_s)/\tan(\pi f/f_s))^{2n})\). See butterworth_lowpass() on precision, and Signal processing: sampled energy, spectra and quantization.

Examples

>>> b, a = butterworth_highpass(1, 250.0, 1000.0)
>>> b.tolist()
[0.5, -0.5]
engineeringkit.signal_processing.butterworth_lowpass(order, cutoff, fs, output='ba')[source]#

Digital Butterworth low-pass filter by the bilinear transform.

Parameters:

  • order (int) – Filter order \(n\), at least 1.

  • cutoff (float) – -3 dB frequency \(f_c\) [Hz], between 0 and \(f_s/2\).

  • fs (float) – Sampling rate [Hz], positive.

  • output ({"ba", "sos"}, optional) – "ba" for numerator and denominator polynomials in \(z^{-1}\); "sos" for a cascade of second-order sections.

Returns:

(b, a) [-], or an (n_sections, 6) array of [b0, b1, b2, 1, a1, a2] rows [-].

Return type:

tuple of ndarray or ndarray

Notes

\(|H(f)|^2 = 1/(1 + (\tan(\pi f/f_s)/\tan(\pi f_c/f_s))^{2n})\): the first \(2n-1\) derivatives vanish at \(f = 0\), and the gain is \(1/\sqrt2\) at \(f_c\) for every order. Above about order 8, or for a cutoff far below \(f_s/2\), the "ba" polynomials lose precision and the filter can come out wrong or unstable; use "sos". See Signal processing: sampled energy, spectra and quantization.

References

S. Butterworth, “On the theory of filter amplifiers,” Wireless Engineer 7 (1930) 536-541; L. B. Jackson, “On the interaction of roundoff noise and dynamic range in digital filters,” Bell Syst. Tech. J. 49 (1970) 159-184.

Examples

>>> b, a = butterworth_lowpass(1, 250.0, 1000.0)
>>> np.round(b, 12).tolist()
[0.5, 0.5]
engineeringkit.signal_processing.chebyshev1_lowpass(order, ripple_db, cutoff, fs, output='ba')[source]#

Digital Chebyshev type I low-pass filter (equiripple passband) by the bilinear transform.

Parameters:

  • order (int) – Filter order \(n\), at least 1.

  • ripple_db (float) – Peak-to-peak passband ripple \(R\) [dB], positive.

  • cutoff (float) – Passband edge [Hz], where the gain last equals \(-R\) dB, between 0 and \(f_s/2\).

  • fs (float) – Sampling rate [Hz], positive.

  • output ({"ba", "sos"}, optional) – Polynomials or second-order sections, as in butterworth_lowpass().

Returns:

(b, a) or second-order sections [-].

Return type:

tuple of ndarray or ndarray

Notes

The analog prototype has \(|H|^2 = 1/(1 + \varepsilon^2 T_n^2(\Omega/\Omega_p))\), with \(T_n\) the Chebyshev polynomial and \(\varepsilon^2 = 10^{R/10} - 1\). \(|T_n| \le 1\) on \([-1, 1]\), so the passband gain swings \(n\) times between 1 and \(1/\sqrt{1+\varepsilon^2}\); outside it \(T_n\) grows as fast as any polynomial of degree \(n\) can, which is why the transition is steeper than Butterworth’s for the same order. Frequencies are mapped by the bilinear transform with the edge prewarped. See Signal processing: sampled energy, spectra and quantization.

References

P. L. Chebyshev, “Théorie des mécanismes connus sous le nom de parallélogrammes,” Mém. Acad. Sci. St.-Pétersbourg, Savants étrangers 7 (1854) 539-568; A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, 3rd ed., Pearson (2010), sec. 7.3.

Examples

>>> b, a = chebyshev1_lowpass(4, 1.0, 100.0, 1000.0)
>>> f, H = frequency_response(b, a, worN=[0.0, 100.0], fs=1000.0)
>>> np.round(20 * np.log10(np.abs(H)), 6).tolist()
[-1.0, -1.0]
engineeringkit.signal_processing.chebyshev2_lowpass(order, attenuation_db, cutoff, fs, output='ba')[source]#

Digital Chebyshev type II (inverse Chebyshev) low-pass filter: flat passband, equiripple stopband.

Parameters:

  • order (int) – Filter order \(n\), at least 1.

  • attenuation_db (float) – Minimum stopband attenuation \(A\) [dB], positive.

  • cutoff (float) – Stopband edge [Hz], where the gain first reaches \(-A\) dB, between 0 and \(f_s/2\).

  • fs (float) – Sampling rate [Hz], positive.

  • output ({"ba", "sos"}, optional) – Polynomials or second-order sections, as in butterworth_lowpass().

Returns:

(b, a) or second-order sections [-].

Return type:

tuple of ndarray or ndarray

Notes

The prototype \(|H|^2 = 1/(1 + 1/(\varepsilon^2 T_n^2(\Omega_s/\Omega)))\) is maximally flat at DC and has zeros on the unit circle in the stopband, where the gain ripples between 0 and \(10^{-A/20}\). See chebyshev1_lowpass() and Signal processing: sampled energy, spectra and quantization.

Examples

>>> b, a = chebyshev2_lowpass(4, 40.0, 150.0, 1000.0)
>>> f, H = frequency_response(b, a, worN=[0.0, 150.0], fs=1000.0)
>>> np.round(20 * np.log10(np.abs(H)), 6).tolist()
[0.0, -40.0]
engineeringkit.signal_processing.circular_convolution(x, h, n=None)[source]#

Circular (periodic) convolution \(y_m = \sum_{k=0}^{N-1} h_k\,x_{(m-k) \bmod N}\), written out.

Parameters:

  • x (array_like) – 1D input samples [signal units].

  • h (array_like) – 1D impulse response [-].

  • n (int, optional) – Period \(N\), at least max(len(x), len(h)); both sequences are zero-padded to it. Defaults to that maximum.

Returns:

The \(N\) output samples [signal units].

Return type:

ndarray

Notes

This is what multiplying two \(N\)-point DFTs computes: \(\mathrm{DFT}(y) = \mathrm{DFT}(x)\,\mathrm{DFT}(h)\). The last \(N_h - 1\) samples of the linear convolution wrap around and add onto the first ones (time aliasing) unless \(N \ge N_x + N_h - 1\), in which case the two agree. See overlap_add() and Signal processing: sampled energy, spectra and quantization.

References

A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, 3rd ed., Pearson (2010), sec. 8.6-8.7.

Examples

>>> circular_convolution([1.0, 2.0, 3.0], [1.0, 1.0]).tolist()
[4.0, 3.0, 5.0]
>>> circular_convolution([1.0, 2.0, 3.0], [1.0, 1.0], n=4).tolist()
[1.0, 3.0, 5.0, 3.0]
engineeringkit.signal_processing.convolution(x, h, mode='full')[source]#

Linear convolution \(y_n = \sum_k h_k x_{n-k}\) (NumPy’s convolve).

Parameters:

  • x (array_like) – Input samples [signal units].

  • h (array_like) – Impulse response [-].

  • mode ({"full", "same", "valid"}, optional) – Output length, as in numpy.convolve().

Returns:

Output samples [signal units].

Return type:

ndarray

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> convolution([1, 2], [1, 1]).tolist()
[1, 3, 2]
engineeringkit.signal_processing.correlation(x, y, mode='full')[source]#

Cross-correlation \(r_k = \sum_n x_{n+k}\,\overline{y_n}\) (SciPy’s correlate).

Parameters:

  • x (array_like) – Samples [signal units].

  • y (array_like) – Samples [signal units].

  • mode ({"full", "same", "valid"}, optional) – Output length, as in scipy.signal.correlate(); "full" returns lags \(-(N_y-1)\) to \(N_x-1\).

Returns:

Correlation at each lag [signal units squared].

Return type:

ndarray

Notes

correlation(x, x) is the unnormalized autocorrelation; divide by \(N\) for the biased estimate whose DFT is the periodogram. See Signal processing: sampled energy, spectra and quantization.

Examples

>>> correlation([1, 2], [1, 2]).tolist()
[2, 5, 2]
engineeringkit.signal_processing.db20(x)[source]#

Amplitude ratio in decibels, \(20\log_{10}|x|\).

Parameters:

x (float or array_like) – Amplitude ratio or gain [-]; complex values use their magnitude.

Returns:

Ratio [dB]; zero gives -inf.

Return type:

float or ndarray

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> float(db20(10.0))
20.0
engineeringkit.signal_processing.dct(x)[source]#

Orthonormal DCT-II \(X_k = \alpha_k \sum_{n=0}^{N-1} x_n \cos\frac{\pi(2n+1)k}{2N}\), written out.

Parameters:

x (array_like) – Real 1D samples [signal units].

Returns:

Coefficients \(X_0..X_{N-1}\) [signal units], with \(\alpha_0 = \sqrt{1/N}\) and \(\alpha_k = \sqrt{2/N}\).

Return type:

ndarray

Notes

The basis is orthonormal, so \(\sum X_k^2 = \sum x_n^2\) and idct() is the transpose. The DCT is the DFT of the record extended by its mirror image, which has no jump at the ends; the DFT’s implicit periodic extension does, so its coefficients fall off far more slowly for a smooth signal. For a first-order Markov process with correlation \(\rho \to 1\) the DCT tends to the Karhunen-Loève transform. Evaluated as an \(N \times N\) matrix product, so \(O(N^2)\); equal to scipy.fft.dct(x, norm="ortho"). See Signal processing: sampled energy, spectra and quantization.

References

N. Ahmed, T. Natarajan and K. R. Rao, “Discrete cosine transform,” IEEE Trans. Comput. C-23 (1974) 90-93, doi:10.1109/T-C.1974.223784.

Examples

>>> (np.round(dct([1.0, 1.0, 1.0, 1.0]), 12) + 0.0).tolist()
[2.0, 0.0, 0.0, 0.0]
engineeringkit.signal_processing.decimate(x, factor, numtaps=None)[source]#

Reduce the sampling rate by an integer factor: anti-alias low-pass FIR filter, then downsample.

Parameters:

  • x (array_like) – 1D samples [signal units].

  • factor (int) – Decimation factor \(M\), at least 1.

  • numtaps (int, optional) – Odd length of the Hamming-window FIR filter, cutoff \(f_s/(2M)\). Defaults to \(20M + 1\).

Returns:

ceil(len(x)/M) samples at rate \(f_s/M\) [signal units], aligned with the input: the filter’s \((N-1)/2\)-sample delay is removed.

Return type:

ndarray

Notes

The filter removes what would alias before downsample() discards samples. Its gain is about -6 dB at the new Nyquist frequency and its stopband about 53 dB down, so keep wanted content below roughly \(0.8 f_s/(2M)\). The record is treated as zero outside its ends, so the first and last \((N-1)/2\) input samples carry a transient. For large \(M\), decimating in stages is cheaper. See Signal processing: sampled energy, spectra and quantization.

References

R. E. Crochiere and L. R. Rabiner, Multirate Digital Signal Processing, Prentice-Hall (1983), ch. 2.

Examples

>>> np.round(decimate(np.ones(200), 4)[20:22], 6).tolist()
[1.0, 1.0]
engineeringkit.signal_processing.dft(x)[source]#

Discrete Fourier transform \(X_k = \sum_{n=0}^{N-1} x_n e^{-j2\pi kn/N}\) (NumPy’s FFT).

Parameters:

x (array_like) – Samples, real or complex [signal units].

Returns:

Unnormalized coefficients \(X_k\) [signal units].

Return type:

ndarray of complex

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> dft([1.0, 1.0, 1.0, 1.0]).real.tolist()
[4.0, 0.0, 0.0, 0.0]
engineeringkit.signal_processing.downsample(x, factor, phase=0)[source]#

Keep every factor-th sample, \(y_m = x_{mM + p}\), with no filtering.

Parameters:

  • x (array_like) – 1D samples [signal units].

  • factor (int) – Downsampling factor \(M\), at least 1.

  • phase (int, optional) – Offset \(p\) of the first kept sample, from 0 to \(M-1\).

Returns:

Samples at rate \(f_s/M\) [signal units].

Return type:

ndarray

Notes

The spectrum of the result is the sum of \(M\) copies of the input spectrum stretched by \(M\) and shifted by multiples of \(2\pi/M\), so any content above \(f_s/(2M)\) aliases. Use decimate() to filter first. See Signal processing: sampled energy, spectra and quantization.

Examples

>>> downsample([0.0, 1.0, 2.0, 3.0, 4.0], 2).tolist()
[0.0, 2.0, 4.0]
engineeringkit.signal_processing.elliptic_lowpass(order, ripple_db, attenuation_db, cutoff, fs, output='ba')[source]#

Digital elliptic (Cauer) low-pass filter: equiripple in both the passband and the stopband.

Parameters:

  • order (int) – Filter order \(n\), at least 1.

  • ripple_db (float) – Peak-to-peak passband ripple \(R\) [dB], positive.

  • attenuation_db (float) – Minimum stopband attenuation \(A\) [dB], above ripple_db.

  • cutoff (float) – Passband edge [Hz], between 0 and \(f_s/2\).

  • fs (float) – Sampling rate [Hz], positive.

  • output ({"ba", "sos"}, optional) – Polynomials or second-order sections, as in butterworth_lowpass().

Returns:

(b, a) or second-order sections [-].

Return type:

tuple of ndarray or ndarray

Notes

The prototype \(|H|^2 = 1/(1 + \varepsilon^2 R_n^2(\Omega/\Omega_p))\) uses Cauer’s elliptic rational function \(R_n\), which ripples between \(\pm 1\) in the passband and stays above \(\sqrt{(10^{A/10}-1)/\varepsilon^2}\) in the stopband. Spreading the error over both bands gives the narrowest transition of any filter of order \(n\) with that ripple and attenuation, at the cost of the most nonlinear phase. See iir_order() and Signal processing: sampled energy, spectra and quantization.

References

W. Cauer, Siebschaltungen, VDI-Verlag, Berlin (1931); W. Cauer, Synthesis of Linear Communication Networks, McGraw-Hill (1958).

Examples

>>> b, a = elliptic_lowpass(4, 0.5, 60.0, 100.0, 1000.0)
>>> f, H = frequency_response(b, a, worN=[100.0], fs=1000.0)
>>> round(float(20 * np.log10(np.abs(H[0]))), 6)
-0.5
engineeringkit.signal_processing.energy(x)[source]#

Discrete signal energy \(\sum_n |x_n|^2\).

Parameters:

x (array_like) – Samples, real or complex [signal units].

Returns:

Energy [signal units squared], without sample-period weighting; multiply by \(1/f_s\) to approximate a continuous integral.

Return type:

float

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> energy([1.0, 2.0])
5.0
engineeringkit.signal_processing.envelope(x)[source]#

Envelope \(|z|\) of a real signal, the magnitude of its analytic signal.

Parameters:

x (array_like) – Real 1D samples [signal units], at least 2.

Returns:

Envelope [signal units].

Return type:

ndarray

Notes

For \(x = a_n\cos(\omega_c n + \phi_n)\) with \(a_n\) and \(\phi_n\) varying slowly compared with the carrier (Bedrosian’s condition), \(|z| = a_n\): an AM detector without a rectifier. See Signal processing: sampled energy, spectra and quantization.

Examples

>>> np.round(envelope(3 * np.cos(2 * np.pi * np.arange(8) / 8)), 12).tolist()
[3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0]
engineeringkit.signal_processing.equiripple_lowpass(numtaps, pass_edge, stop_edge, fs, weight=1.0)[source]#

Parks-McClellan equiripple linear-phase FIR low-pass filter (SciPy’s remez).

Parameters:

  • numtaps (int) – Number of coefficients \(N\), at least 3.

  • pass_edge (float) – Passband edge [Hz].

  • stop_edge (float) – Stopband edge [Hz], with \(0 < f_p < f_s' < f_s/2\).

  • fs (float) – Sampling rate [Hz], positive.

  • weight (float, optional) – Stopband error weight relative to the passband [-], positive; the designed deviations satisfy \(\delta_p = \text{weight}\cdot\delta_s\).

Returns:

Symmetric impulse response [-].

Return type:

ndarray

Notes

Among all symmetric FIR filters of length \(N\), the result minimizes the largest weighted error \(\max |W(f)(H(f) - D(f))|\) over both bands. By the alternation theorem the optimum’s error touches its peak, with alternating sign, at least \((N+1)/2 + 1\) frequencies, and Remez’s exchange algorithm finds them. For the same length it beats any window design, whose largest ripple sits next to the band edge. The transition band is a “don’t care” region. See equiripple_order() and Signal processing: sampled energy, spectra and quantization.

References

T. W. Parks and J. H. McClellan, “Chebyshev approximation for nonrecursive digital filters with linear phase,” IEEE Trans. Circuit Theory 19 (1972) 189-194, doi:10.1109/TCT.1972.1083419.

Examples

>>> h = equiripple_lowpass(41, 100.0, 150.0, 1000.0)
>>> h.size, bool(np.allclose(h, h[::-1]))
(41, True)
engineeringkit.signal_processing.equiripple_order(passband_deviation, stopband_deviation, transition_width, fs)[source]#

Kaiser’s estimate of the Parks-McClellan filter length for given deviations.

Parameters:

  • passband_deviation (float) – Passband deviation \(\delta_p\) [-], in (0, 1).

  • stopband_deviation (float) – Stopband deviation \(\delta_s\) [-], in (0, 1).

  • transition_width (float) – Transition band \(\Delta f\) [Hz], positive and below \(f_s/2\).

  • fs (float) – Sampling rate [Hz], positive.

Returns:

Number of coefficients \(N\).

Return type:

int

Notes

\(N - 1 \approx (-20\log_{10}\sqrt{\delta_p\delta_s} - 13)/(14.6\,\Delta f/f_s)\), rounded up. The estimate is empirical and can be a few taps short; check the design and add taps if needed. Pass weight = passband_deviation / stopband_deviation to equiripple_lowpass(). See Signal processing: sampled energy, spectra and quantization.

References

J. F. Kaiser, “Nonrecursive digital filter design using the I0-sinh window function,” Proc. IEEE Int. Symp. Circuits Syst. (1974) 20-23; Oppenheim and Schafer, Discrete-Time Signal Processing, 3rd ed. (2010), sec. 7.7.

Examples

>>> equiripple_order(0.01, 0.001, 50.0, 1000.0)
52
engineeringkit.signal_processing.fft_frequencies(n, fs)[source]#

Frequency of each DFT bin, \(k f_s/N\), with the upper half negative.

Parameters:

  • n (int) – Transform length \(N\).

  • fs (float) – Sampling rate [Hz], positive.

Returns:

Bin frequencies [Hz] in NumPy’s FFT order.

Return type:

ndarray

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> fft_frequencies(4, 4.0).tolist()
[0.0, 1.0, -2.0, -1.0]
engineeringkit.signal_processing.filter_signal(b, a, x)[source]#

Apply a digital filter, given as polynomials or as second-order sections.

Parameters:

  • b (array_like) – Numerator coefficients in \(z^{-1}\) [-], or, when a is None, an (n_sections, 6) array of second-order sections.

  • a (array_like or None) – Denominator coefficients [-], or None for sections.

  • x (array_like) – Input samples [signal units].

Returns:

Output \(y\), with \(\sum_k a_k y_{n-k} = \sum_k b_k x_{n-k}\) and zero initial state [signal units].

Return type:

ndarray

Notes

The filter runs forward once, so its start-up transient and phase delay remain. Sections are applied one after another (scipy.signal.sosfilt); see butterworth_lowpass() and Signal processing: sampled energy, spectra and quantization.

Examples

>>> filter_signal([1, 1], [1], [1, 0, 0]).tolist()
[1.0, 1.0, 0.0]
engineeringkit.signal_processing.fir_lowpass(numtaps, cutoff, fs, window='hamming')[source]#

Linear-phase FIR low-pass filter by the window method (SciPy’s firwin).

Parameters:

  • numtaps (int) – Number of coefficients \(N\), at least 1.

  • cutoff (float) – Cutoff frequency [Hz], between 0 and \(f_s/2\); the gain there is about one half (-6 dB).

  • fs (float) – Sampling rate [Hz], positive.

  • window (str or tuple, optional) – Window, e.g. "hamming" or ("kaiser", beta); see kaiser_order().

Returns:

Symmetric impulse response [-], scaled to unit gain at DC.

Return type:

ndarray

Notes

The ideal low-pass impulse response \(2(f_c/f_s)\operatorname{sinc}(2 f_c (n - (N-1)/2)/f_s)\) is truncated to \(N\) samples and tapered by the window. Symmetry gives a constant group delay of \((N-1)/2\) samples. See Signal processing: sampled energy, spectra and quantization.

Examples

>>> round(float(fir_lowpass(31, 100.0, 1000.0).sum()), 12)
1.0
engineeringkit.signal_processing.fourier_partial_sum(t, period, a, b, a0=0.0)[source]#

Truncated Fourier series \(a_0 + \sum_{k=1}^{K} [a_k\cos(k\omega t) + b_k\sin(k\omega t)]\).

Parameters:

  • t (float or array_like) – Times [s].

  • period (float) – Period \(2\pi/\omega\) [s], positive.

  • a (array_like) – Cosine and sine coefficients \(a_1..a_K\) and \(b_1..b_K\) [signal units], of equal length.

  • b (array_like) – Cosine and sine coefficients \(a_1..a_K\) and \(b_1..b_K\) [signal units], of equal length.

  • a0 (float, optional) – Mean value [signal units].

Returns:

Partial sum at t [signal units].

Return type:

float or ndarray

Notes

At a jump discontinuity the partial sums overshoot by \(\frac{1}{\pi}\operatorname{Si}(\pi) - \frac12 \approx 0.0895\) of the jump, however large \(K\) is; the overshoot only moves closer to the jump. See Signal processing: sampled energy, spectra and quantization.

References

    1. Gibbs, “Fourier’s series,” Nature 59 (1899) 606.

Examples

>>> float(fourier_partial_sum(0.0, 1.0, [1.0, 0.5], [0.0, 0.0], a0=0.25))
1.75
engineeringkit.signal_processing.frequency_response(b, a=1, worN=512, fs=6.283185307179586)[source]#

Frequency response \(H(e^{j2\pi f/f_s})\) of a digital filter.

Parameters:

  • b (array_like) – Numerator coefficients [-], or second-order sections when a is None.

  • a (array_like or None, optional) – Denominator coefficients [-]; 1 for an FIR filter, None for sections.

  • worN (int or array_like, optional) – Number of frequencies from 0 to \(f_s/2\), or the frequencies themselves [units of fs].

  • fs (float, optional) – Sampling rate [Hz]; the default \(2\pi\) gives frequencies in rad/sample.

Returns:

(frequency, H): frequencies [units of fs] and complex gain [-].

Return type:

tuple of ndarray

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> f, H = frequency_response([0.5, 0.5], 1, worN=[0.0, 250.0], fs=1000.0)
>>> np.round(np.abs(H), 6).tolist()
[1.0, 0.707107]
engineeringkit.signal_processing.goertzel(x, frequency, fs)[source]#

Goertzel’s algorithm: the DFT of x at one frequency from a second-order recursion.

Parameters:

  • x (array_like) – Real samples [signal units].

  • frequency (float or array_like) – Frequency or frequencies to evaluate [Hz]; need not fall on a bin.

  • fs (float) – Sampling rate [Hz], positive.

Returns:

\(X(f) = \sum_n x_n e^{-j2\pi fn/f_s}\) [signal units]; at \(f = kf_s/N\) this is DFT bin \(k\).

Return type:

complex or ndarray of complex

Notes

The recursion \(s_n = x_n + 2\cos\omega\,s_{n-1} - s_{n-2}\), with \(\omega = 2\pi f/f_s\), costs one real multiplication per sample. One extra step with zero input and \(X = e^{-j\omega N}(s_N - e^{-j\omega}s_{N-1})\) finish it. It beats an FFT when fewer than about \(\log_2 N\) bins are wanted, as in a touch-tone (DTMF) decoder. See Signal processing: sampled energy, spectra and quantization.

References

G. Goertzel, “An algorithm for the evaluation of finite trigonometric series,” Amer. Math. Monthly 65 (1958) 34-35.

Examples

>>> X = goertzel([1.0, 0.0, -1.0, 0.0], 1.0, 4.0)
>>> round(float(X.real), 12), round(float(X.imag), 12)
(2.0, 0.0)
engineeringkit.signal_processing.haar_transform(x, levels=None)[source]#

Orthonormal multilevel Haar wavelet transform.

Parameters:

  • x (array_like) – Real 1D samples [signal units]; the length must be divisible by \(2^{levels}\).

  • levels (int, optional) – Number of averaging and differencing stages, at least 1. Defaults to the largest number the length allows.

Returns:

[a_J, d_J, d_{J-1}, ..., d_1]: the coarsest approximation and the details from coarse to fine [signal units], the order used by PyWavelets’ wavedec.

Return type:

list of ndarray

Notes

Each stage maps pairs to \(a = (x_{2k} + x_{2k+1})/\sqrt2\) and \(d = (x_{2k} - x_{2k+1})/\sqrt2\), then repeats on \(a\). The sum of squared coefficients equals that of x. See Signal processing: sampled energy, spectra and quantization.

References

A. Haar, “Zur Theorie der orthogonalen Funktionensysteme,” Math. Ann. 69 (1910) 331-371.

Examples

>>> [c.round(6).tolist() for c in haar_transform([4.0, 2.0, 5.0, 5.0])]
[[8.0], [-2.0], [1.414214, 0.0]]
engineeringkit.signal_processing.hilbert_fir(numtaps, window='hamming')[source]#

Windowed FIR Hilbert transformer, \(h_m = 2/(\pi m)\) for odd \(m\), 0 for even.

Parameters:

  • numtaps (int) – Odd number of coefficients \(N\), at least 3; the taps sit at \(m = n - (N-1)/2\).

  • window (str or tuple, optional) – Window, as in fir_lowpass().

Returns:

Antisymmetric impulse response [-], delay \((N-1)/2\) samples.

Return type:

ndarray

Notes

The ideal \(-j\operatorname{sgn}(\omega)\) response has the impulse response \((1 - \cos\pi m)/(\pi m)\), infinitely long. The truncated, windowed version is a type III linear-phase filter: its gain is zero at DC and Nyquist and close to 1 between. Delay \(x\) by \((N-1)/2\) samples to pair it with the output. See Signal processing: sampled energy, spectra and quantization.

References

A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, 3rd ed., Pearson (2010), sec. 12.4.

Examples

>>> np.round(hilbert_fir(5, window="boxcar") * np.pi, 12).tolist()
[0.0, -2.0, 0.0, 2.0, 0.0]
engineeringkit.signal_processing.hilbert_transform(x)[source]#

Discrete Hilbert transform \(\mathcal{H}x\): every component shifted by -90°.

Parameters:

x (array_like) – Real 1D samples [signal units], at least 2.

Returns:

\(\operatorname{Im} z\) of analytic_signal() [signal units]; \(\cos\omega n\) becomes \(\sin\omega n\).

Return type:

ndarray

Notes

\(\mathcal{H}(\mathcal{H}x) = -x\) for a record with no DC or Nyquist content, and \(x\) and \(\mathcal{H}x\) are orthogonal. See Signal processing: sampled energy, spectra and quantization.

References

D. Hilbert, “Grundzüge einer allgemeinen Theorie der linearen Integralgleichungen (Dritte Mitteilung),” Nachr. Ges. Wiss. Göttingen, Math.-Phys. Kl. (1905) 307-338.

Examples

>>> (np.round(hilbert_transform(np.cos(2 * np.pi * np.arange(4) / 4)), 12) + 0.0).tolist()
[0.0, 1.0, 0.0, -1.0]
engineeringkit.signal_processing.idct(X)[source]#

Inverse of the orthonormal DCT-II (the DCT-III), \(x = C^T X\).

Parameters:

X (array_like) – Real 1D DCT coefficients [signal units].

Returns:

Samples [signal units].

Return type:

ndarray

Notes

Setting the small coefficients to zero before inverting is transform coding, the core of JPEG. See dct() and Signal processing: sampled energy, spectra and quantization.

Examples

>>> np.round(idct([2.0, 0.0, 0.0, 0.0]), 12).tolist()
[1.0, 1.0, 1.0, 1.0]
engineeringkit.signal_processing.idft(X)[source]#

Inverse DFT \(x_n = \frac{1}{N}\sum_k X_k e^{j2\pi kn/N}\) (NumPy’s inverse FFT).

Parameters:

X (array_like) – DFT coefficients [signal units].

Returns:

Samples [signal units].

Return type:

ndarray of complex

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> idft([4.0, 0.0, 0.0, 0.0]).real.tolist()
[1.0, 1.0, 1.0, 1.0]
engineeringkit.signal_processing.iir_order(pass_edge, stop_edge, ripple_db, attenuation_db, fs, family='butterworth')[source]#

Minimum order of a digital low-pass IIR filter that meets a specification.

Parameters:

  • pass_edge (float) – Passband edge \(f_p\) [Hz].

  • stop_edge (float) – Stopband edge \(f_s'\) [Hz], with \(0 < f_p < f_s' < f_s/2\).

  • ripple_db (float) – Largest passband loss \(R\) [dB], positive.

  • attenuation_db (float) – Smallest stopband loss \(A\) [dB], above ripple_db.

  • fs (float) – Sampling rate [Hz], positive.

  • family ({"butterworth", "chebyshev", "elliptic"}, optional) – Approximation; "chebyshev" holds for both types.

Returns:

The smallest order \(n\).

Return type:

int

Notes

With the prewarped selectivity \(k = \tan(\pi f_p/f_s)/\tan(\pi f_s'/f_s)\) and discrimination \(k_1 = \sqrt{(10^{R/10}-1)/(10^{A/10}-1)}\), \(n\) is the ceiling of \(\log k_1/\log k\) (Butterworth), \(\cosh^{-1}(1/k_1)/\cosh^{-1}(1/k)\) (Chebyshev), or \(K(k)K'(k_1)/(K'(k)K(k_1))\) (elliptic), with \(K\) the complete elliptic integral of the first kind and \(K'(k) = K(\sqrt{1-k^2})\). The orders agree with SciPy’s buttord, cheb1ord and ellipord. See Signal processing: sampled energy, spectra and quantization.

References

T. W. Parks and C. S. Burrus, Digital Filter Design, Wiley (1987), ch. 7; A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, 3rd ed., Pearson (2010), sec. 7.3.

Examples

>>> [iir_order(100.0, 150.0, 1.0, 60.0, 1000.0, family) for family in ("butterworth", "chebyshev", "elliptic")]
[17, 9, 6]
engineeringkit.signal_processing.instantaneous_frequency(x, fs)[source]#

Instantaneous frequency \(\frac{1}{2\pi}\frac{d}{dt}\arg z\) from the unwrapped analytic phase.

Parameters:

  • x (array_like) – Real 1D samples [signal units], at least 2.

  • fs (float) – Sampling rate [Hz], positive.

Returns:

len(x) - 1 frequencies [Hz], each between two samples.

Return type:

ndarray

Notes

The phase difference of consecutive samples, times \(f_s/2\pi\). It is meaningful for a single component; two tones give a frequency that oscillates between them. See Signal processing: sampled energy, spectra and quantization.

Examples

>>> x = np.cos(2 * np.pi * 50 * np.arange(100) / 1000)
>>> np.round(instantaneous_frequency(x, 1000.0)[40:42], 6).tolist()
[50.0, 50.0]
engineeringkit.signal_processing.interpolate(x, factor, numtaps=None)[source]#

Raise the sampling rate by an integer factor: upsample, then low-pass filter with gain \(L\).

Parameters:

  • x (array_like) – 1D samples [signal units].

  • factor (int) – Interpolation factor \(L\), at least 1.

  • numtaps (int, optional) – Odd length of the Hamming-window FIR filter, cutoff \(f_s/2\) of the input. Defaults to \(20L + 1\).

Returns:

len(x) * L samples at rate \(L f_s\) [signal units], aligned with the input.

Return type:

ndarray

Notes

The filter is the windowed ideal interpolator \(h_n = w_n \operatorname{sinc}(n/L)\), which is 1 at \(n = 0\) and 0 at every other multiple of \(L\), so the original samples pass through unchanged and only the new ones are computed. It is Shannon’s sinc reconstruction, truncated to \(N\) taps; the images left by upsample() are attenuated by about 53 dB. See Signal processing: sampled energy, spectra and quantization.

References

R. W. Schafer and L. R. Rabiner, “A digital signal processing approach to interpolation,” Proc. IEEE 61 (1973) 692-702, doi:10.1109/PROC.1973.9150.

Examples

>>> y = interpolate([0.0, 1.0, 0.0, -1.0], 2)
>>> np.round(y[::2], 12).tolist()
[0.0, 1.0, 0.0, -1.0]
engineeringkit.signal_processing.inverse_haar_transform(coefficients)[source]#

Rebuild a signal from haar_transform() coefficients.

Parameters:

coefficients (sequence of array_like) – [a_J, d_J, ..., d_1], each detail twice as long as the one before it, the first detail as long as a_J [signal units].

Returns:

Samples [signal units].

Return type:

ndarray

Notes

Inverts each stage with \(x_{2k} = (a + d)/\sqrt2\) and \(x_{2k+1} = (a - d)/\sqrt2\). Zeroing small details before inverting is the simplest wavelet compression or denoising. See Signal processing: sampled energy, spectra and quantization.

Examples

>>> inverse_haar_transform([[8.0], [-2.0], [np.sqrt(2), 0.0]]).round(12).tolist()
[4.0, 2.0, 5.0, 5.0]
engineeringkit.signal_processing.kaiser_order(attenuation_db, transition_width, fs)[source]#

Kaiser’s formulas for the FIR length and window shape that meet a stopband attenuation.

Parameters:

  • attenuation_db (float) – Required stopband attenuation (and passband ripple) \(A = -20\log_{10}\delta\) [dB], positive.

  • transition_width (float) – Width of the transition band \(\Delta f\) [Hz], positive and below \(f_s/2\).

  • fs (float) – Sampling rate [Hz], positive.

Returns:

(numtaps, beta): number of coefficients (int) and Kaiser window parameter \(\beta\) [-], for window=("kaiser", beta).

Return type:

tuple

Notes

All parameters are scalar-only: the result sizes one filter.

With \(\Delta\omega = 2\pi\Delta f/f_s\), \(N - 1 = \lceil (A - 7.95)/(2.285\,\Delta\omega) \rceil\) and \(\beta = 0.1102(A - 8.7)\) for \(A > 50\), \(0.5842(A-21)^{0.4} + 0.07886(A-21)\) for \(21 \le A \le 50\), and 0 below. The formulas are empirical and usually land within about a decibel of \(A\). See Signal processing: sampled energy, spectra and quantization.

References

J. F. Kaiser, “Digital filters,” ch. 7 in F. F. Kuo and J. F. Kaiser (eds.), System Analysis by Digital Computer, Wiley (1966); J. F. Kaiser, “Nonrecursive digital filter design using the I0-sinh window function,” Proc. IEEE Int. Symp. Circuits Syst. (1974) 20-23.

Examples

>>> numtaps, beta = kaiser_order(60.0, 50.0, 1000.0)
>>> numtaps, round(beta, 4)
(74, 5.6533)
engineeringkit.signal_processing.levinson_durbin(r, order)[source]#

Levinson-Durbin recursion: AR coefficients of every order up to order from autocorrelations.

Parameters:

  • r (array_like) – Real autocorrelation \(r_0, r_1, \ldots\) [signal units squared], with \(r_0 > 0\) and at least order + 1 values.

  • order (int) – Model order \(p\), at least 1.

Returns:

Coefficients \([1, a_1..a_p]\), noise variance \(E_p\), reflection coefficients and prediction errors.

Return type:

ARModel

Notes

Starting from \(E_0 = r_0\), order \(m\) adds the reflection coefficient \(k_m = -(r_m + \sum_{i=1}^{m-1} a_i r_{m-i})/E_{m-1}\), updates \(a_i \leftarrow a_i + k_m a_{m-i}\) with \(a_m = k_m\), and shrinks the error to \(E_m = (1 - k_m^2)E_{m-1}\). The result solves the Yule-Walker equations \(\sum_k a_k r_{|i-k|} = -r_i\), \(i = 1..p\). A \(|k_m| \ge 1\) means the sequence is not a valid (positive definite) autocorrelation. See Signal processing: sampled energy, spectra and quantization.

References

N. Levinson, “The Wiener RMS (root mean square) error criterion in filter design and prediction,” J. Math. Phys. 25 (1947) 261-278; J. Durbin, “The fitting of time-series models,” Rev. Int. Stat. Inst. 28 (1960) 233-244, doi:10.2307/1401322.

Examples

>>> model = levinson_durbin([1.0, 0.5, 0.25], 2)
>>> (np.round(model.coefficients, 12) + 0.0).tolist(), round(model.noise_variance, 12)
([1.0, -0.5, 0.0], 0.75)
engineeringkit.signal_processing.lms_filter(x, d, n_taps, mu, initial_weights=None)[source]#

Widrow-Hoff least-mean-squares adaptive FIR filter.

Parameters:

  • x (array_like) – Real input samples [signal units].

  • d (array_like) – Desired output, same length [signal units].

  • n_taps (int) – Filter length \(L\), between 1 and len(x).

  • mu (float) – Step size [1/signal units squared], positive.

  • initial_weights (array_like, optional) – Starting weights, length n_taps [-]; zeros by default.

Returns:

Output, error, final weights and the weight history.

Return type:

LMSResult

Notes

For each sample, with \(u_n = [x_n, \ldots, x_{n-L+1}]\): \(y_n = w^T u_n\), \(e_n = d_n - y_n\) and \(w \leftarrow w + \mu e_n u_n\). The mean weights converge to the Wiener solution when \(0 < \mu < 2/\lambda_{max}(R)\); in practice \(\mu \ll 2/(L\,E[x^2])\) keeps the excess error small. See Signal processing: sampled energy, spectra and quantization.

References

B. Widrow and M. E. Hoff, “Adaptive switching circuits,” IRE WESCON Convention Record, part 4 (1960) 96-104.

Examples

>>> rng = np.random.default_rng(0)
>>> x = rng.normal(size=2000)
>>> result = lms_filter(x, 0.5 * x, 1, 0.05)
>>> round(float(result.weights[0]), 6)
0.5
engineeringkit.signal_processing.moving_average(x, n)[source]#

Moving average over n samples, \(y_m = \frac{1}{n}\sum_{k=0}^{n-1} x_{m+k}\).

Parameters:

  • x (array_like) – 1D input samples [signal units].

  • n (int) – Window length, between 1 and len(x).

Returns:

The len(x) - n + 1 fully overlapped averages [signal units].

Return type:

ndarray

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> moving_average([1.0, 2.0, 3.0], 2).tolist()
[1.5, 2.5]
engineeringkit.signal_processing.nyquist_frequency(fs)[source]#

Nyquist frequency \(f_s/2\), the highest frequency the samples represent unambiguously.

Parameters:

fs (float or array_like) – Sampling rate [Hz], positive.

Returns:

Nyquist frequency [Hz].

Return type:

float or ndarray

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> float(nyquist_frequency(100))
50.0
engineeringkit.signal_processing.overlap_add(x, h, block_length=None)[source]#

Linear convolution of a long record with an FIR filter by FFT blocks, overlap-add.

Parameters:

  • x (array_like) – 1D input samples [signal units].

  • h (array_like) – 1D impulse response of length \(M\) [-].

  • block_length (int, optional) – Input samples per block \(L\), at least 1. Defaults to the largest \(L\) for which \(L + M - 1\) fits in the smallest power of two of at least \(4M\).

Returns:

The \(N_x + M - 1\) samples of convolution(x, h) [signal units], real when both inputs are real.

Return type:

ndarray

Notes

The record is cut into blocks of \(L\) samples. Each block and \(h\) are zero-padded to \(N_{FFT} \ge L + M - 1\), so that the circular convolution of the block equals its linear one, multiplied in the frequency domain, and the \(L + M - 1\) output samples are added into place, overlapping the next block’s by \(M - 1\). The cost per output sample grows as \(\log N_{FFT}\) rather than \(M\). See circular_convolution() and Signal processing: sampled energy, spectra and quantization.

References

T. G. Stockham Jr., “High-speed convolution and correlation,” AFIPS Conf. Proc. 28 (1966 Spring Joint Computer Conf.) 229-233, doi:10.1145/1464182.1464209.

Examples

>>> np.round(overlap_add([1.0, 2.0, 3.0, 4.0], [1.0, 1.0], block_length=2), 12).tolist()
[1.0, 3.0, 5.0, 7.0, 4.0]
engineeringkit.signal_processing.power(x)[source]#

Mean signal power \(\overline{|x|^2}\).

Parameters:

x (array_like) – Samples, real or complex [signal units].

Returns:

Power [signal units squared].

Return type:

float

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> power([1.0, -1.0])
1.0
engineeringkit.signal_processing.power_spectral_density(x, fs)[source]#

One-sided periodogram \(P_k = \frac{2|X_k|^2}{f_s N}\) (SciPy’s periodogram).

Parameters:

  • x (array_like) – Samples [signal units].

  • fs (float) – Sampling rate [Hz], positive.

Returns:

(frequency, psd): bin frequencies [Hz] and power spectral density [signal units squared per Hz]; its sum times the bin spacing is the signal variance.

Return type:

tuple of ndarray

Notes

The periodogram’s standard deviation equals its mean at every bin, so a longer record sharpens it without smoothing it; see welch_psd(). SciPy removes the mean before transforming. See Signal processing: sampled energy, spectra and quantization.

References

A. Schuster, “On the investigation of hidden periodicities with application to a supposed 26 day period of meteorological phenomena,” Terr. Magn. 3 (1898) 13-41.

engineeringkit.signal_processing.prewarp_frequency(frequency, fs)[source]#

Analog frequency that the bilinear transform maps onto a given digital frequency.

Parameters:

  • frequency (float or array_like) – Digital frequency \(f\) [Hz], from 0 to below \(f_s/2\).

  • fs (float) – Sampling rate [Hz], positive.

Returns:

Analog frequency \(f_a = (f_s/\pi)\tan(\pi f/f_s)\) [Hz].

Return type:

float or ndarray

Notes

The bilinear transform with \(K = 2f_s\) takes the analog axis \(s = j2\pi f_a\) onto the unit circle at \(f = (f_s/\pi)\arctan(\pi f_a/f_s)\), compressing all of \([0, \infty)\) into \([0, f_s/2)\). Designing the analog prototype at \(f_a\) puts its edge at \(f\) after the transform. \(f_a \approx f\) only for \(f \ll f_s\). See bilinear_transform() and Signal processing: sampled energy, spectra and quantization.

Examples

>>> round(float(prewarp_frequency(250.0, 1000.0)), 6)
318.309886
engineeringkit.signal_processing.quantization_step(vmin, vmax, bits)[source]#

Step \(q = (v_{max} - v_{min})/2^N\) of a uniform mid-rise ADC.

Parameters:

  • vmin (float) – Input range [V or any amplitude unit], vmax > vmin.

  • vmax (float) – Input range [V or any amplitude unit], vmax > vmin.

  • bits (int) – Resolution \(N\), between 1 and 52.

Returns:

Quantization step [same unit as vmin].

Return type:

float

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> quantization_step(-1.0, 1.0, 3)
0.25
engineeringkit.signal_processing.quantize(x, vmin, vmax, bits)[source]#

Ideal uniform mid-rise quantization, saturating outside the ADC range.

Parameters:

  • x (float or array_like) – Input samples [V or any amplitude unit].

  • vmin (float) – Input range [same unit], vmax > vmin.

  • vmax (float) – Input range [same unit], vmax > vmin.

  • bits (int) – Resolution \(N\), between 1 and 52.

Returns:

One of the \(2^N\) reconstruction levels at the bin centres \(v_{min} + (k + 1/2)q\) [same unit].

Return type:

float or ndarray

Notes

In range, the error never exceeds \(q/2\); for a busy signal it is close to uniform on \([-q/2, q/2]\), with power \(q^2/12\). See Signal processing: sampled energy, spectra and quantization.

References

W. R. Bennett, “Spectra of quantized signals,” Bell Syst. Tech. J. 27 (1948) 446-472.

Examples

>>> quantize([-0.9, 0.1], -1.0, 1.0, 3).tolist()
[-0.875, 0.125]
engineeringkit.signal_processing.radix2_fft(x, return_multiplications=False)[source]#

Recursive decimation-in-time radix-2 FFT, written out to show the butterflies.

Parameters:

  • x (array_like) – Samples, real or complex [signal units]; the length must be a power of two.

  • return_multiplications (bool, optional) – Also return the number of complex twiddle multiplications.

Returns:

DFT coefficients, equal to dft() up to rounding, and, if requested, the multiplication count \((N/2)\log_2 N\).

Return type:

ndarray of complex or tuple

Notes

Each stage combines the DFTs \(E_k\) and \(O_k\) of the even and odd samples as \(X_k = E_k + W^k O_k\) and \(X_{k+N/2} = E_k - W^k O_k\), with \(W = e^{-j2\pi/N}\). The count includes the trivial twiddle \(W^0 = 1\). See Signal processing: sampled energy, spectra and quantization.

References

J. W. Cooley and J. W. Tukey, “An algorithm for the machine calculation of complex Fourier series,” Math. Comp. 19 (1965) 297-301.

Examples

>>> X, count = radix2_fft([1.0, 0.0, 0.0, 0.0], return_multiplications=True)
>>> X.real.tolist(), count
([1.0, 1.0, 1.0, 1.0], 4)
engineeringkit.signal_processing.real_cepstrum(x, n_fft=None)[source]#

Real cepstrum \(c_n = \mathrm{IDFT}(\ln|\mathrm{DFT}(x)|)\), indexed by quefrency \(n/f_s\).

Parameters:

  • x (array_like) – Real 1D samples [signal units].

  • n_fft (int, optional) – Transform length, at least len(x); zero-padding reduces the wrap-around of the cepstrum. Defaults to len(x).

Returns:

Real, even cepstrum of length n_fft [-]; sample \(n\) is quefrency \(n/f_s\) [s].

Return type:

ndarray

Notes

An echo \(x_n + \alpha x_{n-d}\) multiplies the spectrum by \(1 + \alpha e^{-j\omega d}\), a ripple in \(\ln|X|\) with “frequency” \(d\); its log expands as \(\sum_m (-1)^{m+1}(\alpha^m/m)\cos(m\omega d)\), so for \(|\alpha| < 1\) the cepstrum adds \(\alpha/2\) at \(n = d\) (and smaller peaks at multiples of \(d\)) to that of \(x\). Convolution becomes addition, which is how speech is split into pitch and vocal tract. A zero in the spectrum has no logarithm and raises ValueError. See Signal processing: sampled energy, spectra and quantization.

References

B. P. Bogert, M. J. R. Healy and J. W. Tukey, “The quefrency alanysis of time series for echoes: cepstrum, pseudo-autocovariance, cross-cepstrum and saphe cracking,” in M. Rosenblatt (ed.), Proceedings of the Symposium on Time Series Analysis, Wiley (1963) 209-243.

Examples

>>> x = np.zeros(64)
>>> x[0], x[10] = 1.0, 0.5
>>> round(float(real_cepstrum(x)[10]), 6)
0.25
engineeringkit.signal_processing.rms(x)[source]#

Root-mean-square value \(\sqrt{\overline{|x|^2}}\).

Parameters:

x (array_like) – Samples, real or complex [signal units].

Returns:

RMS value [signal units].

Return type:

float

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> float(rms([3.0, -3.0]))
3.0
engineeringkit.signal_processing.sampling_period(fs)[source]#

Sampling period \(T = 1/f_s\).

Parameters:

fs (float or array_like) – Sampling rate [Hz], positive.

Returns:

Sampling period [s].

Return type:

float or ndarray

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> float(sampling_period(10))
0.1
engineeringkit.signal_processing.sinc_reconstruct(samples, fs, t)[source]#

Whittaker-Shannon interpolation \(x(t) = \sum_n x_n \operatorname{sinc}(f_s t - n)\).

Parameters:

  • samples (array_like) – Samples \(x_n\) taken at \(t_n = n/f_s\) [signal units].

  • fs (float) – Sampling rate [Hz], positive.

  • t (array_like) – Times at which to evaluate the reconstruction [s].

Returns:

Reconstructed signal at t [signal units].

Return type:

ndarray

Notes

The sum is truncated to the samples supplied, so it is exact at the sample instants but only approximate between them near the ends of the record. See Signal processing: sampled energy, spectra and quantization.

References

C. E. Shannon, “Communication in the presence of noise,” Proc. IRE 37 (1949) 10-21.

Examples

>>> sinc_reconstruct([0.0, 1.0, 0.0], 10.0, [0.1]).tolist()
[1.0]
engineeringkit.signal_processing.snr_db(signal_values, noise)[source]#

Signal-to-noise ratio \(10\log_{10}(P_s/P_n)\) from sample arrays.

Parameters:

  • signal_values (array_like) – Signal samples, real or complex [signal units].

  • noise (array_like) – Noise samples of the same shape [signal units].

Returns:

SNR [dB]; inf for zero noise and -inf for zero signal.

Return type:

float

Raises:

ValueError – If the shapes differ or both powers are zero.

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> float(snr_db([1.0, -1.0], [0.1, -0.1]))
20.0
engineeringkit.signal_processing.stft(x, fs, window_length, hop=None, window='hann')[source]#

Short-time Fourier transform \(X_m(f_k) = \sum_n x_{n+mH}\,w_n\,e^{-j2\pi kn/L}\).

Parameters:

  • x (array_like) – Real samples [signal units].

  • fs (float) – Sampling rate [Hz], positive.

  • window_length (int) – Frame length \(L\) [samples], between 2 and len(x).

  • hop (int, optional) – Frame step \(H\) [samples], at least 1. Defaults to \(L/4\) (rounded down, at least 1).

  • window (str or tuple, optional) – Analysis window, as in windowed_spectrum(); Gabor’s choice is ("gaussian", std).

Returns:

Frame centre times [s], one-sided bin frequencies [Hz] and the complex coefficients [signal units].

Return type:

STFTResult

Notes

A window of duration \(\Delta t\) resolves frequencies only to about \(1/\Delta t\). For RMS widths of \(|w|^2\) in time and of \(|W|^2\) in frequency, \(\sigma_t\sigma_f \ge 1/(4\pi)\), with equality only for a Gaussian. Frames are not padded, so the first and last centres lie \(L/2\) samples inside the record. See Signal processing: sampled energy, spectra and quantization.

References

D. Gabor, “Theory of communication,” J. IEE 93, part III (1946) 429-457.

Examples

>>> result = stft(np.ones(16), 16.0, 8, hop=4, window="boxcar")
>>> result.summary()["frames"], float(abs(result.values[0, 0]))
(3, 8.0)
engineeringkit.signal_processing.upsample(x, factor)[source]#

Insert \(L - 1\) zeros after every sample, \(y_{nL} = x_n\).

Parameters:

  • x (array_like) – 1D samples [signal units].

  • factor (int) – Upsampling factor \(L\), at least 1.

Returns:

len(x) * L samples at rate \(L f_s\) [signal units].

Return type:

ndarray

Notes

\(Y(e^{j\omega}) = X(e^{j\omega L})\): the spectrum is compressed by \(L\) and \(L - 1\) images of it appear up to the new Nyquist frequency. interpolate() removes them. See Signal processing: sampled energy, spectra and quantization.

Examples

>>> upsample([1.0, 2.0], 3).tolist()
[1.0, 0.0, 0.0, 2.0, 0.0, 0.0]
engineeringkit.signal_processing.welch_psd(x, fs, segment_length=256, overlap=None, window='hann')[source]#

Welch’s averaged, modified periodogram, written out segment by segment.

Parameters:

  • x (array_like) – Real samples [signal units].

  • fs (float) – Sampling rate [Hz], positive.

  • segment_length (int, optional) – Samples per segment \(L\), between 2 and len(x).

  • overlap (int, optional) – Samples shared by consecutive segments, from 0 to \(L-1\). Defaults to \(L/2\) (rounded down).

  • window (str or tuple, optional) – Window applied to each segment, as in windowed_spectrum().

Returns:

(frequency, psd): one-sided bin frequencies [Hz] and power spectral density [signal units squared per Hz].

Return type:

tuple of ndarray

Notes

Each segment’s periodogram is \(|\mathrm{FFT}(w\,x)|^2/(f_s\sum w^2)\), doubled for every bin but DC and Nyquist, and the segments are averaged. Averaging \(K\) nearly independent segments divides the variance by about \(K\) and widens the resolution to about \(f_s/L\). No detrending is applied, so the result matches scipy.signal.welch(..., detrend=False). See Signal processing: sampled energy, spectra and quantization.

References

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

Examples

>>> f, psd = welch_psd(np.ones(8), 8.0, segment_length=4, window="boxcar")
>>> psd.tolist()
[0.5, 0.0, 0.0]
engineeringkit.signal_processing.wiener_filter(x, d, n_taps)[source]#

FIR Wiener filter: the weights that minimize \(E[(d_n - \sum_k w_k x_{n-k})^2]\).

Parameters:

  • x (array_like) – Real input samples [signal units].

  • d (array_like) – Desired output, same length [signal units].

  • n_taps (int) – Filter length \(L\), between 1 and len(x).

Returns:

Weights \(w_0..w_{L-1}\) [-], to be applied as filter_signal(w, 1, x).

Return type:

ndarray

Notes

Solves \(Rw = p\) with the biased sample estimates \(r_k = \frac1N\sum_n x_n x_{n-k}\) and \(p_k = \frac1N\sum_n d_n x_{n-k}\), by Levinson recursion (scipy.linalg.solve_toeplitz). The minimum mean squared error is \(E[d^2] - p^T w\). Signals are assumed stationary. See Signal processing: sampled energy, spectra and quantization.

References

N. Wiener, Extrapolation, Interpolation, and Smoothing of Stationary Time Series, MIT Press and Wiley (1949); N. Levinson, “The Wiener RMS (root mean square) error criterion in filter design and prediction,” J. Math. Phys. 25 (1947) 261-278.

Examples

>>> x = np.array([1.0, -1.0, 2.0, 0.5, -1.5, 1.0])
>>> np.round(wiener_filter(x, 3 * x, 2), 12).tolist()
[3.0, 0.0]
engineeringkit.signal_processing.windowed_spectrum(x, fs, window='hann', n_fft=None)[source]#

One-sided magnitude \(|\mathrm{FFT}(w\,x)|\) of a windowed record.

Parameters:

  • x (array_like) – Real samples [signal units].

  • fs (float) – Sampling rate [Hz], positive.

  • window (str or tuple, optional) – Window name accepted by scipy.signal.get_window(), e.g. "boxcar", "hann", "hamming", "blackman" or ("kaiser", beta).

  • n_fft (int, optional) – Transform length, at least len(x); the windowed record is zero-padded to interpolate the spectrum. Defaults to len(x).

Returns:

(frequency, magnitude): bin frequencies [Hz] and unnormalized magnitude [signal units].

Return type:

tuple of ndarray

Notes

A tone between bins leaks into its neighbours through the window’s sidelobes. The window does not change the bin spacing \(f_s/N\); zero-padding only samples the same spectrum more finely. See Signal processing: sampled energy, spectra and quantization.

References

R. B. Blackman and J. W. Tukey, “The measurement of power spectra from the point of view of communications engineering,” Bell Syst. Tech. J. 37 (1958) 185-282 and 485-569.

Examples

>>> f, magnitude = windowed_spectrum([1.0, 1.0, 1.0, 1.0], 4.0, window="boxcar")
>>> magnitude.tolist()
[4.0, 0.0, 0.0]
engineeringkit.signal_processing.yule_walker(x, order, demean=True)[source]#

Yule-Walker estimate of an AR model from a record, solved by levinson_durbin().

Parameters:

  • x (array_like) – Real 1D samples [signal units], longer than order.

  • order (int) – Model order \(p\), at least 1.

  • demean (bool, optional) – Subtract the record’s mean first.

Returns:

The fitted model; model.psd(fs) is the AR spectrum.

Return type:

ARModel

Notes

Uses the biased autocorrelation \(r_k = \frac1N\sum_n x_n x_{n+k}\), which is always positive definite, so the model is always stable. Its spectrum is smooth and can resolve two close peaks from a record too short for the periodogram, but a wrong order gives spurious or missing peaks. See Signal processing: sampled energy, spectra and quantization.

References

G. U. Yule, “On a method of investigating periodicities in disturbed series, with special reference to Wolfer’s sunspot numbers,” Phil. Trans. R. Soc. A 226 (1927) 267-298, doi:10.1098/rsta.1927.0007; G. Walker, “On periodicity in series of related terms,” Proc. R. Soc. A 131 (1931) 518-532.

Examples

>>> rng = np.random.default_rng(0)
>>> e = rng.normal(size=20000)
>>> x = np.zeros_like(e)
>>> for n in range(1, e.size):
...     x[n] = 0.9 * x[n - 1] + e[n]
>>> round(float(yule_walker(x, 1).coefficients[1]), 2)
-0.9

Utilities#

Supporting numerics for engineeringkit.signal_processing’s systems/ modules.

engineeringkit.signal_processing.utils.db20(x)[source]#

Amplitude ratio in decibels, \(20\log_{10}|x|\).

Parameters:

x (float or array_like) – Amplitude ratio or gain [-]; complex values use their magnitude.

Returns:

Ratio [dB]; zero gives -inf.

Return type:

float or ndarray

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> float(db20(10.0))
20.0
engineeringkit.signal_processing.utils.energy(x)[source]#

Discrete signal energy \(\sum_n |x_n|^2\).

Parameters:

x (array_like) – Samples, real or complex [signal units].

Returns:

Energy [signal units squared], without sample-period weighting; multiply by \(1/f_s\) to approximate a continuous integral.

Return type:

float

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> energy([1.0, 2.0])
5.0
engineeringkit.signal_processing.utils.power(x)[source]#

Mean signal power \(\overline{|x|^2}\).

Parameters:

x (array_like) – Samples, real or complex [signal units].

Returns:

Power [signal units squared].

Return type:

float

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> power([1.0, -1.0])
1.0
engineeringkit.signal_processing.utils.rms(x)[source]#

Root-mean-square value \(\sqrt{\overline{|x|^2}}\).

Parameters:

x (array_like) – Samples, real or complex [signal units].

Returns:

RMS value [signal units].

Return type:

float

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> float(rms([3.0, -3.0]))
3.0
engineeringkit.signal_processing.utils.snr_db(signal_values, noise)[source]#

Signal-to-noise ratio \(10\log_{10}(P_s/P_n)\) from sample arrays.

Parameters:

  • signal_values (array_like) – Signal samples, real or complex [signal units].

  • noise (array_like) – Noise samples of the same shape [signal units].

Returns:

SNR [dB]; inf for zero noise and -inf for zero signal.

Return type:

float

Raises:

ValueError – If the shapes differ or both powers are zero.

Notes

Inputs are raw SI numbers with no dimensional metadata; see Signal processing: sampled energy, spectra and quantization.

Examples

>>> float(snr_db([1.0, -1.0], [0.1, -0.1]))
20.0

Visualizers#

Matplotlib visualizers for engineeringkit.signal_processing.

engineeringkit.signal_processing.visualizers.plot_magnitude_response(b, a=1, fs=6.283185307179586, n=2048, ax=None, label=None)[source]#

Plot a digital filter’s gain in decibels from DC to the Nyquist frequency.

Parameters:

  • b (array_like) – Filter polynomials, or second-order sections in b with a=None, as in frequency_response().

  • a (array_like) – Filter polynomials, or second-order sections in b with a=None, as in frequency_response().

  • fs (float, optional) – Sampling rate [Hz]; the default \(2\pi\) gives rad/sample.

  • n (int, optional) – Number of frequencies.

  • ax (matplotlib.axes.Axes, optional) – Axes to draw on; a new figure is created if omitted.

  • label (str, optional) – Legend label.

Return type:

matplotlib.axes.Axes

engineeringkit.signal_processing.visualizers.plot_pole_zero(b, a, ax=None)[source]#

Plot the zeros (o) and poles (x) of \(B(z)/A(z)\) with the unit circle.

Parameters:

  • b (array_like) – Numerator and denominator coefficients in \(z^{-1}\) [-].

  • a (array_like) – Numerator and denominator coefficients in \(z^{-1}\) [-].

  • ax (matplotlib.axes.Axes, optional) – Axes to draw on; a new figure is created if omitted.

Return type:

matplotlib.axes.Axes

On this page
Show Source