Source code for PyDynamic.misc.filterstuff

r"""This module is a collection of functions which are related to filter design

This module contains the following functions:

* :func:`db`: Calculation of decibel values :math:`20\log_{10}(x)` for a vector of
  values
* :func:`ua`: Shortcut for calculation of unwrapped angle of complex values
* :func:`grpdelay`: Calculation of the group delay of a digital filter
* :func:`mapinside`: Maps the roots of polynomial with coefficients :math:`a`
  to the unit circle
* :func:`kaiser_lowpass`: Design of a FIR lowpass filter using the window technique
  with a Kaiser window.
* :func:`isstable`: Determine whether an IIR filter with certain coefficients is stable
* :func:`savitzky_golay`: Smooth (and optionally differentiate) data with a
  Savitzky-Golay filter
"""

__all__ = [
    "db",
    "ua",
    "grpdelay",
    "mapinside",
    "kaiser_lowpass",
    "isstable",
    "savitzky_golay",
]

import numpy as np


[docs] def db(vals): r"""Calculation of decibel values :math:`20\log_{10}(x)` for a vector of values""" return 20 * np.log10(np.abs(vals))
[docs] def ua(vals): """Shortcut for calculation of unwrapped angle of complex values""" return np.unwrap(np.angle(vals))
[docs] def grpdelay(b, a, Fs, nfft=512): """Calculation of the group delay of a digital filter Parameters ---------- b: ndarray IIR filter numerator coefficients a: ndarray IIR filter denominator coefficients Fs: float sampling frequency of the filter nfft: int number of FFT bins Returns ------- group_delay: np.ndarray group delay values frequencies: ndarray frequencies at which the group delay is calculated References ---------- * Smith, online book [Smith]_ """ Na = len(a) - 1 Nb = len(b) - 1 c = np.convolve(b, a[::-1]) # c(z) = b(z)*a(1/z)*z^(-oa) cr = c * np.arange(Na + Nb + 1) # derivative of c wrt 1/z num = np.fft.fft(cr, 2 * nfft) den = np.fft.fft(c, 2 * nfft) tol = 1e-12 polebins = np.nonzero(abs(den) < tol) for p in polebins: num[p] = 0.0 den[p] = 1.0 gd = np.real(num / den) - Na f = np.arange(0.0, 2 * nfft - 1) / (2 * nfft) * Fs f = f[: nfft + 1] gd = gd[: len(f)] return gd, f
[docs] def mapinside(a): """Maps the roots of polynomial to the unit circle. Maps the roots of polynomial with coefficients :math:`a` to the unit circle. Parameters ---------- a: ndarray polynomial coefficients Returns ------- a: ndarray polynomial coefficients with all roots inside or on the unit circle """ v = np.roots(a) inds = np.nonzero(abs(v) > 1) v[inds] = 1 / np.conj(v[inds]) return np.poly(v)
[docs] def kaiser_lowpass(L, fcut, Fs, beta=8.0): """Design of a FIR lowpass filter using the window technique with a Kaiser window. This method uses a Kaiser window. Filters of that type are often used as FIR low-pass filters due to their linear phase. Parameters ---------- L: int filter order (window length) fcut: float desired cut-off frequency Fs: float sampling frequency beta: float scaling parameter for the Kaiser window Returns ------- blow: ndarray FIR filter coefficients shift: int delay of the filter (in samples) """ from scipy.signal import firwin if np.mod(L, 2) == 0: L += 1 blow = firwin(L, 2 * fcut / Fs, window=("kaiser", beta)) shift = L / 2 return blow, shift
[docs] def isstable(b, a, ftype="digital"): """Determine whether an IIR filter with certain coefficients is stable Determine whether IIR filter with coefficients `b` and `a` is stable by checking the roots of the polynomial `a`. Parameters ---------- b : ndarray Filter numerator coefficients. These are only part of the input parameters for compatibility reasons (especially with MATLAB code). During the computation they are actually not used. a : ndarray Filter denominator coefficients. ftype : string, optional Filter type. 'digital' if in discrete-time (default) and 'analog' if in continuous-time. Returns ------- stable : bool Whether filter is stable or not. """ v = np.roots(a) # Check if all of the roots of the polynomial made from the filter coefficients... if ftype == "digital": # ... lie inside the unit circle in discrete time. return not np.any(np.abs(v) > 1.0) elif ftype == "analog": # ... are non-negative in continuous time. return not np.any(np.real(v) < 0) else: raise ValueError( f"isstable: filter type 'ftype={ftype}'unknown. Please " f"check documentation and choose a valid assignment for " f"'ftype' or leave out to fallback to default." )
[docs] def savitzky_golay(y, window_size, order, deriv=0, delta=1.0): """Smooth (and optionally differentiate) data with a Savitzky-Golay filter The Savitzky-Golay filter removes high frequency noise from data. It has the advantage of preserving the original shape and features of the signal better than other types of filtering approaches, such as moving averages techniques. Source obtained from scipy cookbook (online), downloaded 2013-09-13 Parameters ---------- y: ndarray, shape (N,) the values of the time history of the signal window_size: int the length of the window. Must be an odd integer number order: int the order of the polynomial used in the filtering. Must be less then `window_size` - 1. deriv: int, optional The order of the derivative to compute. This must be a nonnegative integer. The default is 0, which means to filter the data without differentiating. delta: float, optional The spacing of the samples to which the filter will be applied. This is only used if deriv > 0. This includes a factor :math:`n! / h^n`, where :math:`n` is represented by `deriv` and :math:`1/h` by `delta`. Returns ------- ys: ndarray, shape (N,) the smoothed signal (or it's n-th derivative). Notes ----- The Savitzky-Golay is a type of low-pass filter, particularly suited for smoothing noisy data. The main idea behind this approach is to make for each point a least-square fit with a polynomial of high order over a odd-sized window centered at the point. References ---------- * Savitzky et al. [Savitzky]_ * Numerical Recipes [NumRec]_ """ from math import factorial try: window_size = np.abs(np.int(window_size)) order = np.abs(np.int(order)) except ValueError: raise ValueError("window_size and order have to be of type int") if window_size % 2 != 1 or window_size < 1: raise TypeError("window_size size must be a positive odd number") if window_size < order + 2: raise TypeError("window_size is too small for the polynomials order") order_range = range(order + 1) half_window = (window_size - 1) // 2 # precompute coefficients b = np.mat( [[k**i for i in order_range] for k in range(-half_window, half_window + 1)] ) m = np.linalg.pinv(b).A[deriv] * factorial(deriv) / delta**deriv # pad the signal at the extremes with # values taken from the signal itself firstvals = y[0] - np.abs(y[1 : half_window + 1][::-1] - y[0]) lastvals = y[-1] + np.abs(y[-half_window - 1 : -1][::-1] - y[-1]) y = np.concatenate((firstvals, y, lastvals)) return np.convolve(m[::-1], y, mode="valid")