Source code for PyDynamic.model_estimation.fit_filter

# -*- coding: utf-8 -*-
"""
The module :mod:`PyDynamic.model_estimation.fit_filter` contains several functions to
carry out a least-squares fit to a given complex frequency response and the design of
digital deconvolution filters by least-squares fitting to the reciprocal of a given
frequency response each with associated uncertainties.

This module contains the following functions:

* :func:`LSIIR`: Least-squares IIR filter fit to a given frequency response
* :func:`LSFIR`: Least-squares fit of a digital FIR filter to a given frequency
  response
* :func:`invLSFIR`: Least-squares fit of a digital FIR filter to the reciprocal of a
  given frequency response.
* :func:`invLSFIR_unc`: Design of FIR filter as fit to reciprocal of frequency response
  values with uncertainty
* :func:`invLSFIR_uncMC`: Design of FIR filter as fit to reciprocal of frequency
  response values with uncertainty via Monte Carlo
* :func:`invLSIIR`: Design of a stable IIR filter as fit to reciprocal of frequency
  response values
* :func:`invLSIIR_unc`: Design of a stable IIR filter as fit to reciprocal of frequency
  response values with uncertainty

"""
import warnings

import numpy as np
import scipy.signal as dsp

from ..misc.filterstuff import grpdelay, mapinside

__all__ = [
    "LSIIR",
    "LSFIR",
    "invLSFIR",
    "invLSFIR_unc",
    "invLSIIR",
    "invLSIIR_unc",
    "invLSFIR_uncMC"
    ]


def _fitIIR(
    Hvals: np.ndarray,
    tau: int,
    w: np.ndarray,
    E: np.ndarray,
    Na: int,
    Nb: int,
    inv: bool = False
    ):
    """The actual fitting routing for the least-squares IIR filter.

    Parameters
    ----------
        Hvals :  (M,) np.ndarray
            (complex) frequency response values
        tau : integer
            initial estimate of time delay
        w : np.ndarray
            :math:`2 * np.pi * f / Fs`
        E : np.ndarray
            :math:`np.exp(-1j * np.dot(w[:, np.newaxis], Ns.T))`
        Nb : int
            numerator polynomial order
        Na : int
            denominator polynomial order
        inv : bool, optional
            If True the least-squares fitting is performed for the reciprocal, if False
            for the actual frequency response

    Returns
    -------
        b, a : IIR filter coefficients as numpy arrays
    """
    exponent = -1 if inv else 1
    Ea = E[:, 1:Na + 1]
    Eb = E[:, :Nb + 1]
    Htau = np.exp(-1j * w * tau) * Hvals ** exponent
    HEa = np.dot(np.diag(Htau), Ea)
    D = np.hstack((HEa, -Eb))
    Tmp1 = np.real(np.dot(np.conj(D.T), D))
    Tmp2 = np.real(np.dot(np.conj(D.T), -Htau))
    ab = np.linalg.lstsq(Tmp1, Tmp2, rcond=None)[0]
    a = np.hstack((1.0, ab[:Na]))
    b = ab[Na:]
    return b, a


[docs]def LSIIR(Hvals, Nb, Na, f, Fs, tau=0, justFit=False): """Least-squares IIR filter fit to a given frequency response. This method uses Gauss-Newton non-linear optimization and pole mapping for filter stabilization Parameters ---------- Hvals: numpy array of (complex) frequency response values of shape (M,) Nb: integer numerator polynomial order Na: integer denominator polynomial order f: numpy array of frequencies at which Hvals is given of shape (M,) Fs: sampling frequency tau: integer initial estimate of time delay justFit: boolean, when true then no stabilization is carried out Returns ------- b,a: IIR filter coefficients as numpy arrays tau: filter time delay in samples References ---------- * Eichstädt et al. 2010 [Eichst2010]_ * Vuerinckx et al. 1996 [Vuer1996]_ """ print("\nLeast-squares fit of an order %d digital IIR filter" % max(Nb, Na)) print("to a frequency response given by %d values.\n" % len(Hvals)) w = 2 * np.pi * f / Fs Ns = np.arange(0, max(Nb, Na) + 1)[:, np.newaxis] E = np.exp(-1j * np.dot(w[:, np.newaxis], Ns.T)) b, a = _fitIIR(Hvals, tau, w, E, Na, Nb, inv=False) if justFit: print("Calculation done. No stabilization requested.") if np.count_nonzero(np.abs(np.roots(a)) > 1) > 0: print("Obtained filter is NOT stable.") sos = np.sum( np.abs((dsp.freqz(b, a, 2 * np.pi * f / Fs)[1] - Hvals) ** 2)) print("Final sum of squares = %e" % sos) tau = 0 return b, a, tau if np.count_nonzero(np.abs(np.roots(a)) > 1) > 0: stable = False else: stable = True maxiter = 50 astab = mapinside(a) run = 1 while stable is not True and run < maxiter: g1 = grpdelay(b, a, Fs)[0] g2 = grpdelay(b, astab, Fs)[0] tau = np.ceil(tau + np.median(g2 - g1)) b, a = _fitIIR(Hvals, tau, w, E, Na, Nb, inv=False) if np.count_nonzero(np.abs(np.roots(a)) > 1) > 0: astab = mapinside(a) else: stable = True run = run + 1 if np.count_nonzero(np.abs(np.roots(a)) > 1) > 0: print("Caution: The algorithm did NOT result in a stable IIR filter!") print( "Maybe try again with a higher value of tau0 or a higher filter " "order?") print( "Least squares fit finished after %d iterations (tau=%d)." % (run, tau)) Hd = dsp.freqz(b, a, 2 * np.pi * f / Fs)[1] Hd = Hd * np.exp(1j * 2 * np.pi * f / Fs * tau) res = np.hstack( (np.real(Hd) - np.real(Hvals), np.imag(Hd) - np.imag(Hvals))) rms = np.sqrt(np.sum(res ** 2) / len(f)) print("Final rms error = %e \n\n" % rms) return b, a, int(tau)
[docs]def LSFIR(H, N, tau, f, Fs, Wt=None): """ Least-squares fit of a digital FIR filter to a given frequency response. Parameters ---------- H : (complex) frequency response values of shape (M,) N : FIR filter order tau : delay of filter f : frequencies of shape (M,) Fs : sampling frequency of digital filter Wt : (optional) vector of weights of shape (M,) or shape (M,M) Returns ------- filter coefficients bFIR (ndarray) of shape (N+1,) """ print("\nLeast-squares fit of an order %d digital FIR filter to the" % N) print("reciprocal of a frequency response given by %d values.\n" % len(H)) H = H[:, np.newaxis] w = 2 * np.pi * f / Fs w = w[:, np.newaxis] ords = np.arange(N + 1)[:, np.newaxis] ords = ords.T E = np.exp(-1j * np.dot(w, ords)) if Wt is not None: if len(np.shape(Wt)) == 2: # is matrix weights = np.diag(Wt) else: weights = np.eye(len(f)) * Wt X = np.vstack( [np.real(np.dot(weights, E)), np.imag(np.dot(weights, E))]) else: X = np.vstack([np.real(E), np.imag(E)]) H = H * np.exp(1j * w * tau) iRI = np.vstack([np.real(1.0 / H), np.imag(1.0 / H)]) bFIR, res = np.linalg.lstsq(X, iRI)[:2] if not isinstance(res, np.ndarray): print( "Calculation of FIR filter coefficients finished with residual " "norm %e" % res) return np.reshape(bFIR, (N + 1,))
[docs]def invLSFIR(H, N, tau, f, Fs, Wt=None): """ Least-squares fit of a digital FIR filter to the reciprocal of a given frequency response. Parameters ---------- H: np.ndarray of shape (M,) and dtype complex frequency response values N: int FIR filter order tau: float delay of filter f: np.ndarray of shape (M,) frequencies Fs: float sampling frequency of digital filter Wt: np.ndarray of shape (M,) - optional vector of weights Returns ------- bFIR: np.ndarray of shape (N,) filter coefficients References ---------- * Elster and Link [Elster2008]_ .. see_also ::mod::`PyDynamic.uncertainty.propagate_filter.FIRuncFilter` """ print("\nLeast-squares fit of an order %d digital FIR filter to the" % N) print("reciprocal of a frequency response given by %d values.\n" % len(H)) H = H[:, np.newaxis] # extend to matrix-like for simplified algebra w = 2 * np.pi * f / Fs # set up radial frequencies w = w[:, np.newaxis] ords = np.arange(N + 1)[:, np.newaxis] # set up design matrix ords = ords.T E = np.exp(-1j * np.dot(w, ords)) if Wt is not None: # set up weighted design matrix if necessary if len(np.shape(Wt)) == 2: # is matrix weights = np.diag(Wt) else: weights = np.eye(len(f)) * Wt X = np.vstack([np.real(np.dot(weights, E)), np.imag(np.dot(weights, E))]) else: X = np.vstack([np.real(E), np.imag(E)]) Hs = H * np.exp(1j * w * tau) # apply time delay for improved fit quality iRI = np.vstack([np.real(1.0 / Hs), np.imag(1.0 / Hs)]) bFIR, res = np.linalg.lstsq(X, iRI)[:2] # the actual fitting if (not isinstance(res, np.ndarray)) or (len(res) == 1): # summarise results print( "Calculation of FIR filter coefficients finished with residual " "norm %e" % res ) Hd = dsp.freqz(bFIR, 1, 2 * np.pi * f / Fs)[1] Hd = Hd * np.exp(1j * 2 * np.pi * f / Fs * tau) res = np.hstack((np.real(Hd) - np.real(H), np.imag(Hd) - np.imag(H))) rms = np.sqrt(np.sum(res ** 2) / len(f)) print("Final rms error = %e \n\n" % rms) return bFIR.flatten()
[docs]def invLSFIR_unc(H, UH, N, tau, f, Fs, wt=None, verbose=True, trunc_svd_tol=None): """Design of FIR filter as fit to reciprocal of frequency response values with uncertainty Least-squares fit of a digital FIR filter to the reciprocal of a frequency response for which associated uncertainties are given for its real and imaginary part. Uncertainties are propagated using a truncated svd and linear matrix propagation. Parameters ---------- H: np.ndarray of shape (M,) frequency response values UH: np.ndarray of shape (2M,2M) uncertainties associated with the real and imaginary part N: int FIR filter order tau: float delay of filter f: np.ndarray of shape (M,) frequencies Fs: float sampling frequency of digital filter wt: np.ndarray of shape (2M,) - optional array of weights for a weighted least-squares method verbose: bool, optional whether to print statements to the command line trunc_svd_tol: float lower bound for singular values to be considered for pseudo-inverse Returns ------- b: np.ndarray of shape (N+1,) filter coefficients of shape Ub: np.ndarray of shape (N+1,N+1) uncertainties associated with b References ---------- * Elster and Link [Elster2008]_ """ if verbose: print("\nLeast-squares fit of an order %d digital FIR filter to the" % N) print("reciprocal of a frequency response given by %d values" % len(H)) print("and propagation of associated uncertainties.") # Step 1: Propagation of uncertainties to reciprocal of frequency response runs = 10000 Nf = len(f) if not len(H) == UH.shape[0]: # Assume that H is given as complex valued frequency response. RI = np.hstack((np.real(H), np.imag(H))) else: RI = H.copy() H = H[:Nf] + 1j * H[Nf:] HRI = np.random.multivariate_normal(RI, UH, runs) # random draws of real,imag of # freq response values omtau = 2 * np.pi * f / Fs * tau # Vectorized Monte Carlo for propagation to inverse absHMC = HRI[:, :Nf] ** 2 + HRI[:, Nf:] ** 2 HiMC = np.hstack(((HRI[:, :Nf] * np.tile(np.cos(omtau), (runs, 1)) + HRI[:, Nf:] * np.tile(np.sin(omtau), (runs, 1))) / absHMC, (HRI[:, Nf:] * np.tile(np.cos(omtau), (runs, 1)) - HRI[:, :Nf] * np.tile(np.sin(omtau), (runs, 1))) / absHMC)) UiH = np.cov(HiMC, rowvar=False) # Step 2: Fit filter coefficients and evaluate uncertainties if isinstance(wt, np.ndarray): if wt.shape != np.diag(UiH).shape[0]: raise ValueError("User-defined weighting has wrong dimension.") else: wt = np.ones(2 * Nf) E = np.exp(-1j * 2 * np.pi * np.dot(f[:, np.newaxis] / Fs, np.arange(N + 1)[:, np.newaxis].T)) X = np.vstack((np.real(E), np.imag(E))) X = np.dot(np.diag(wt), X) Hm = H * np.exp(1j * 2 * np.pi * f / Fs * tau) Hri = np.hstack((np.real(1.0 / Hm), np.imag(1.0 / Hm))) u, s, v = np.linalg.svd(X, full_matrices=False) if isinstance(trunc_svd_tol, float): s[s < trunc_svd_tol] = 0.0 StSInv = np.zeros_like(s) StSInv[s > 0] = s[s > 0] ** (-2) M = np.dot(np.dot(np.dot(v.T, np.diag(StSInv)), np.diag(s)), u.T) bFIR = np.dot(M, Hri[:, np.newaxis]) # actual fitting UbFIR = np.dot(np.dot(M, UiH), M.T) # evaluation of uncertainties bFIR = bFIR.flatten() if verbose: Hd = dsp.freqz(bFIR, 1, 2 * np.pi * f / Fs)[1] Hd = Hd * np.exp(1j * 2 * np.pi * f / Fs * tau) res = np.hstack((np.real(Hd) - np.real(H), np.imag(Hd) - np.imag(H))) rms = np.sqrt(np.sum(res ** 2) / len(f)) print("Final rms error = %e \n\n" % rms) return bFIR, UbFIR
[docs]def invLSFIR_uncMC(H, UH, N, tau, f, Fs, verbose=True): """Design of FIR filter as fit to reciprocal of frequency response values with uncertainty Least-squares fit of a FIR filter to the reciprocal of a frequency response for which associated uncertainties are given for its real and imaginary parts. Uncertainties are propagated using a Monte Carlo method. This method may help in cases where the weighting matrix or the Jacobian are ill-conditioned, resulting in false uncertainties associated with the filter coefficients. Parameters ---------- H: np.ndarray of shape (M,) and dtype complex frequency response values UH: np.ndarray of shape (2M,2M) uncertainties associated with the real and imaginary part of H N: int FIR filter order tau: int time delay of filter in samples f: np.ndarray of shape (M,) frequencies corresponding to H Fs: float sampling frequency of digital filter verbose: bool, optional whether to print statements to the command line Returns ------- b: np.ndarray of shape (N+1,) filter coefficients of shape Ub: np.ndarray of shape (N+1, N+1) uncertainties associated with b References ---------- * Elster and Link [Elster2008]_ """ if verbose: print("\nLeast-squares fit of an order %d digital FIR filter to the" % N) print("reciprocal of a frequency response given by %d values" % len(H)) print("and propagation of associated uncertainties.") # Step 1: Propagation of uncertainties to reciprocal of frequency response runs = 10000 HRI = np.random.multivariate_normal(np.hstack((np.real(H), np.imag(H))), UH, runs) # Step 2: Fitting the filter coefficients E = np.exp(-1j * 2 * np.pi * np.dot(f[:, np.newaxis] / Fs, np.arange(N + 1)[:, np.newaxis].T)) X = np.vstack((np.real(E), np.imag(E))) Nf = len(f) bF = np.zeros((N + 1, runs)) resn = np.zeros((runs,)) for k in range(runs): Hk = HRI[k, :Nf] + 1j * HRI[k, Nf:] Hkt = Hk * np.exp(1j * 2 * np.pi * f / Fs * tau) iRI = np.hstack([np.real(1.0 / Hkt), np.imag(1.0 / Hkt)]) bF[:, k], res = np.linalg.lstsq(X, iRI)[:2] resn[k] = np.linalg.norm(res) bFIR = np.mean(bF, axis=1) UbFIR = np.cov(bF, rowvar=True) return bFIR, UbFIR
[docs]def invLSIIR(Hvals, Nb, Na, f, Fs, tau, justFit=False, verbose=True): """Design of a stable IIR filter as fit to reciprocal of frequency response values Least-squares fit of a digital IIR filter to the reciprocal of a given set of frequency response values using the equation-error method and stabilization by pole mapping and introduction of a time delay. Parameters ---------- Hvals: np.ndarray of shape (M,) and dtype complex frequency response values. Nb: int order of IIR numerator polynomial. Na: int order of IIR denominator polynomial. f: np.ndarray of shape (M,) frequencies corresponding to Hvals Fs: float sampling frequency for digital IIR filter. tau: float initial estimate of time delay for filter stabilization. justFit: bool if True then no stabilization is carried out. verbose: bool If True print some more detail about input parameters. Returns ------- b : np.ndarray The IIR filter numerator coefficient vector in a 1-D sequence a : np.ndarray The IIR filter denominator coefficient vector in a 1-D sequence tau : int time delay (in samples) References ---------- * Eichstädt, Elster, Esward, Hessling [Eichst2010]_ """ from numpy import count_nonzero, roots, ceil, median if verbose: print( "\nLeast-squares fit of an order %d digital IIR filter to the" % max(Nb, Na) ) print("reciprocal of a frequency response given by %d values.\n" % len(Hvals)) w = 2 * np.pi * f / Fs Ns = np.arange(0, max(Nb, Na) + 1)[:, np.newaxis] E = np.exp(-1j * np.dot(w[:, np.newaxis], Ns.T)) bi, ai = _fitIIR(Hvals, tau, w, E, Na, Nb, inv=True) if justFit: # no uncertainty evaluation return bi, ai if count_nonzero(abs(roots(ai)) > 1) > 0: stable = False else: stable = True maxiter = 50 astab = mapinside(ai) # stabilise filter run = 1 while stable is not True and run < maxiter: # shift delay such that filter # becomes stable g1 = grpdelay(bi, ai, Fs)[0] g2 = grpdelay(bi, astab, Fs)[0] tau = ceil(tau + median(g2 - g1)) bi, ai = _fitIIR(Hvals, tau, w, E, Na, Nb, inv=True) if count_nonzero(abs(roots(ai)) > 1) > 0: astab = mapinside(ai) else: stable = True run = run + 1 if count_nonzero(abs(roots(ai)) > 1) > 0 and verbose: print("Caution: The algorithm did NOT result in a stable IIR filter!") print( "Maybe try again with a higher value of tau0 or a higher filter order?" ) if verbose: print("Least squares fit finished after %d iterations (tau=%d).\n" % (run, tau)) Hd = dsp.freqz(bi, ai, 2 * np.pi * f / Fs)[1] Hd = Hd * np.exp(1j * 2 * np.pi * f / Fs * tau) res = np.hstack((np.real(Hd) - np.real(Hvals), np.imag(Hd) - np.imag(Hvals))) rms = np.sqrt(np.sum(res ** 2) / len(f)) print("Final rms error = %e \n\n" % rms) return bi, ai, int(tau)
[docs]def invLSIIR_unc(H, UH, Nb, Na, f, Fs, tau=0): """Design of stabel IIR filter as fit to reciprocal of given frequency response with uncertainty Least-squares fit of a digital IIR filter to the reciprocal of a given set of frequency response values with given associated uncertainty. Propagation of uncertainties is carried out using the Monte Carlo method. Parameters ---------- H: np.ndarray of shape (M,) and dtype complex frequency response values. UH: np.ndarray of shape (2M,2M) uncertainties associated with real and imaginary part of H Nb: int order of IIR numerator polynomial. Na: int order of IIR denominator polynomial. f: np.ndarray of shape (M,) frequencies corresponding to H Fs: float sampling frequency for digital IIR filter. tau: float initial estimate of time delay for filter stabilization. Returns ------- b,a: np.ndarray IIR filter coefficients tau: int time delay (in samples) Uba: np.ndarray of shape (Nb+Na+1, Nb+Na+1) uncertainties associated with [a[1:],b] References ---------- * Eichstädt, Elster, Esward and Hessling [Eichst2010]_ .. seealso:: :mod:`PyDynamic.uncertainty.propagate_filter.IIRuncFilter` :mod:`PyDynamic.model_estimation.fit_filter.invLSIIR` """ runs = 1000 print("\nLeast-squares fit of an order %d digital IIR filter to the" % max(Nb, Na)) print("reciprocal of a frequency response given by %d values.\n" % len(H)) print( "Uncertainties of the filter coefficients are evaluated using\n" "the GUM S2 Monte Carlo method with %d runs.\n" % runs ) # Step 1: Propagation of uncertainties to frequency response HRI = np.random.multivariate_normal(np.hstack((np.real(H), np.imag(H))), UH, runs) HH = HRI[:, : len(f)] + 1j * HRI[:, len(f) :] # Step 2: Fit filter and evaluate uncertainties (Monte Carlo method) AB = np.zeros((runs, Nb + Na + 1)) Tau = np.zeros((runs,)) for k in range(runs): bi, ai, Tau[k] = invLSIIR(HH[k, :], Nb, Na, f, Fs, tau, verbose=False) AB[k, :] = np.hstack((ai[1:], bi)) bi = np.mean(AB[:, Na:], axis=0) ai = np.hstack((np.array([1.0]), np.mean(AB[:, :Na], axis=0))) Uab = np.cov(AB, rowvar=False) tau = np.mean(Tau) return bi, ai, tau, Uab