Source code for PyDynamic.identification.fit_filter

# -*- coding: utf-8 -*-
"""
This module contains several functions to carry out a least-squares fit to a
given complex frequency response.

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


.. deprecated:: 1.2.71
          The module *identification* will be combined with the module
          *deconvolution* and renamed to *model_estimation* in the next
          major release 2.0.0. From then on you should only use the new module
          *model_estimation* instead.
"""
import warnings

import numpy as np
import scipy.signal as dsp

from ..misc.filterstuff import grpdelay, mapinside

__all__ = ['LSIIR', 'LSFIR']

warnings.simplefilter('default')
warnings.warn(
    "The module *identification* will be combined with the module "
    "*deconvolution* and renamed to *model_estimation* in the "
    "next major release 2.0.0. From then on you should only use "
    "the new module *model_estimation* instead.",
    DeprecationWarning)


def fitIIR(Hvals, tau, w, E, Na, Nb):
    """The actual fitting routing for the least-squares IIR filter.

    Parameters
    ----------
        Hvals:   numpy array of (complex) frequency response values of shape (M,)
        tau:     integer initial estimate of time delay
        w:       numpy array :math:`2 * np.pi * f / Fs`
        E:       numpy array :math:`np.exp(-1j * np.dot(w[:, np.newaxis],
                 Ns.T))`
        Nb:      integer numerator polynomial order
        Na:      integer denominator polynomial order

    Returns
    -------
        b, a:    IIR filter coefficients as numpy arrays
    """
    Ea = E[:, 1:Na + 1]
    Eb = E[:, :Nb + 1]
    Htau = np.exp(-1j * w * tau) * Hvals
    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) 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) 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,))