Source code for PyDynamic.identification.fit_filter

# -*- coding: utf-8 -*-
'''
.. moduleauthor:: Sascha Eichstaedt (sascha.eichstaedt@ptb.de)
'''

from ..misc.filterstuff import grpdelay, mapinside

import numpy as np
import scipy.signal as dsp


[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 frequency response values Nb: integer numerator polynomial order Na: integer denominator polynomial order f: numpy array of frequencies at which Hvals is given 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)) def fitIIR(Hvals,tau,E,Na,Nb): # The actual fitting routine 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)[0] a = np.hstack((1.0,ab[:Na])) b = ab[Na:] return b,a b,a = fitIIR(Hvals,tau,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!=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,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. :param H: frequency response values :param N: FIR filter order :param tau: delay of filter :param f: frequencies :param Fs: sampling frequency of digital filter :param Wt: (optional) vector of weights :type H: ndarray :type N: int :type tau: int :type f: ndarray :type Fs: float :returns: filter coefficients bFIR (ndarray) of shape (N,) """ 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 not Wt == 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,))