Skip to content
Snippets Groups Projects
Commit 2a3a874a authored by Eric Kooistra's avatar Eric Kooistra
Browse files

Original file from https://github.com/chipmuenk.

parent 16d4a0d4
No related branches found
No related tags found
1 merge request!374Add plot_iir_filter_analysis() based on LTF-IIR-allgemein.ipynbcode from...
Show whitespace changes
Inline Side-by-side
applications/lofar2/model/pfb_os/dsp_fpga_lib.py 0 → 100644
+ 1199
− 0
View file @ 2a3a874a
# -*- coding: utf-8 -*-
"""
Created on Mon Apr 30 10:29:42 2012
dsp_fpga_lib (Version 8)
@author: Muenker_2
from: https://github.com/chipmuenk/dsp
"""
import sys
import string # needed for remezord?
import numpy as np
import numpy.ma as ma
from numpy import pi, asarray, absolute, sqrt, log10, arctan, ceil, hstack, mod
import scipy.signal as sig
from scipy import __version__ as sci_version
from scipy import special # needed for remezord
import scipy.spatial.distance as sc_dist
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib import patches
__version__ = "0.9"
def versions():
print("Python version:", ".".join(map(str, sys.version_info[:3])))
print("Numpy:", np.__version__)
print("Scipy:", sci_version)
print("Matplotlib:", mpl.__version__, mpl.get_backend())
mpl_rc = {'lines.linewidth' : 1.5,
'lines.markersize' : 8, # markersize in points
'text.color' : 'black',
'font.family' : 'sans-serif',#'serif',
'font.style' : 'normal',
#'mathtext.fontset' : 'stixsans',#'stix',
#'mathtext.fallback_to_cm' : True,
'mathtext.default' : 'it',
'font.size' : 12,
'legend.fontsize' : 12,
'axes.labelsize' : 12,
'axes.titlesize' : 14,
'axes.linewidth' : 1, # linewidth for coordinate system
'axes.formatter.use_mathtext': True, # use mathtext for scientific notation.
#
'axes.grid' : True,
'grid.color' : '#222222',
'axes.facecolor' : 'white',
'axes.labelcolor' : 'black',
'axes.edgecolor' : 'black',
'grid.linestyle' : ':',
'grid.linewidth' : 0.5,
#
'xtick.top' : False, # mpl >= 2.0
'xtick.direction' : 'out',
'ytick.direction' : 'out',
'xtick.color' : 'black',
'ytick.color' : 'black',
#
'figure.figsize' : (7,4), # default figure size in inches
'figure.facecolor' : 'white',
'figure.edgecolor' : '#808080',
'figure.dpi' : 100,
'savefig.dpi' : 100,
'savefig.facecolor' : 'white',
'savefig.edgecolor' : 'white',
'savefig.bbox' : 'tight',
'savefig.pad_inches' : 0,
'hatch.color' : '#808080', # mpl >= 2.0
'hatch.linewidth' : 0.5, # mpl >= 2.0
'animation.html' : 'jshtml' # javascript, mpl >= 2.1
}
mpl_rc_33 = {'mathtext.fallback' : 'cm'} # new since mpl 3.3
mpl_rc_32 = {'mathtext.fallback_to_cm': True} # deprecated since mpl 3.3
if mpl.__version__ < "3.3":
mpl_rc.update(mpl_rc_32) # lower than matplotlib 3.3
else:
mpl_rc.update(mpl_rc_33)
plt.rcParams.update(mpl_rc) # define plot properties
def H_mag(zaehler, nenner, z, lim):
""" Calculate magnitude of H(z) or H(s) in polynomial form at the complex
coordinate z = x, 1j * y (skalar or array)
The result is clipped at lim."""
# limvec = lim * np.ones(len(z))
try: len(zaehler)
except TypeError:
z_val = abs(zaehler) # zaehler is a scalar
else:
z_val = abs(np.polyval(zaehler,z)) # evaluate zaehler at z
try: len(nenner)
except TypeError:
n_val = nenner # nenner is a scalar
else:
n_val = abs(np.polyval(nenner,z))
return np.minimum((z_val/n_val),lim)
#----------------------------------------------
# from scipy.sig.signaltools.py:
def cmplx_sort(p):
"sort roots based on magnitude."
p = np.asarray(p)
if np.iscomplexobj(p):
indx = np.argsort(abs(p))
else:
indx = np.argsort(p)
return np.take(p, indx, 0), indx
# adapted from scipy.signal.signaltools.py:
# TODO: comparison of real values has several problems (5 * tol ???)
def unique_roots(p, tol=1e-3, magsort = False, rtype='min', rdist='euclidian'):
"""
Determine unique roots and their multiplicities from a list of roots.
Parameters
----------
p : array_like
The list of roots.
tol : float, default tol = 1e-3
The tolerance for two roots to be considered equal. Default is 1e-3.
magsort: Boolean, default = False
When magsort = True, use the root magnitude as a sorting criterium (as in
the version used in numpy < 1.8.2). This yields false results for roots
with similar magniudes (e.g. on the unit circle) but is signficantly
faster for a large number of roots (factor 20 for 500 double roots.)
rtype : {'max', 'min, 'avg'}, optional
How to determine the returned root if multiple roots are within
`tol` of each other.
- 'max' or 'maximum': pick the maximum of those roots (magnitude ?).
- 'min' or 'minimum': pick the minimum of those roots (magnitude?).
- 'avg' or 'mean' : take the average of those roots.
- 'median' : take the median of those roots
dist : {'manhattan', 'euclid'}, optional
How to measure the distance between roots: 'euclid' is the euclidian
distance. 'manhattan' is less common, giving the
sum of the differences of real and imaginary parts.
Returns
-------
pout : list
The list of unique roots, sorted from low to high (only for real roots).
mult : list
The multiplicity of each root.
Notes
-----
This utility function is not specific to roots but can be used for any
sequence of values for which uniqueness and multiplicity has to be
determined. For a more general routine, see `numpy.unique`.
Examples
--------
>>> vals = [0, 1.3, 1.31, 2.8, 1.25, 2.2, 10.3]
>>> uniq, mult = sp.signal.unique_roots(vals, tol=2e-2, rtype='avg')
Check which roots have multiplicity larger than 1:
>>> uniq[mult > 1]
array([ 1.305])
Find multiples of complex roots on the unit circle:
>>> vals = np.roots(1,2,3,2,1)
uniq, mult = sp.signal.unique_roots(vals, rtype='avg')
"""
def manhattan(a,b):
"""
Manhattan distance between a and b
"""
return ma.abs(a.real - b.real) + ma.abs(a.imag - b.imag)
def euclid(a,b):
"""
Euclidian distance between a and b
"""
return ma.abs(a - b)
if rtype in ['max', 'maximum']:
comproot = ma.max # nanmax ignores nan's
elif rtype in ['min', 'minimum']:
comproot = ma.min # nanmin ignores nan's
elif rtype in ['avg', 'mean']:
comproot = ma.mean # nanmean ignores nan's
# elif rtype == 'median':
else:
raise TypeError(rtype)
if rdist in ['euclid', 'euclidian']:
dist_roots = euclid
elif rdist in ['rect', 'manhattan']:
dist_roots = manhattan
else:
raise TypeError(rdist)
mult = [] # initialize list for multiplicities
pout = [] # initialize list for reduced output list of roots
p = np.atleast_1d(p) # convert p to at least 1D array
tol = abs(tol)
if len(p) == 0: # empty argument, return empty lists
return pout, mult
elif len(p) == 1: # scalar input, return arg with multiplicity = 1
pout = p
mult = [1]
return pout, mult
else:
sameroots = [] # temporary list for roots within the tolerance
pout = p[np.isnan(p)].tolist() # copy nan elements to pout as list
mult = len(pout) * [1] # generate a list with a "1" for each nan
#p = ma.masked_array(p[~np.isnan(p)]) # delete nan elements, convert to ma
p = np.ma.masked_where(np.isnan(p), p) # only masks nans, preferrable?
if np.iscomplexobj(p) and not magsort:
for i in range(len(p)): # p[i] is current root under test
if not p[i] is ma.masked: # has current root been "deleted" yet?
tolarr = dist_roots(p[i], p[i:]) < tol # test against itself and
# subsequent roots, giving a multiplicity of at least one
mult.append(np.sum(tolarr)) # multiplicity = number of "hits"
sameroots = p[i:][tolarr] # pick the roots within the tolerance
p[i:] = ma.masked_where(tolarr, p[i:]) # and "delete" (mask) them
pout.append(comproot(sameroots)) # avg/mean/max of mult. root
else:
p,indx = cmplx_sort(p)
indx = -1
curp = p[0] + 5 * tol # needed to avoid "self-detection" ?
for k in range(len(p)):
tr = p[k]
# if dist_roots(tr, curp) < tol:
if abs(tr - curp) < tol:
sameroots.append(tr)
curp = comproot(sameroots) # not correct for 'avg'
# of multiple (N > 2) root !
pout[indx] = curp
mult[indx] += 1
else:
pout.append(tr)
curp = tr
sameroots = [tr]
indx += 1
mult.append(1)
return np.array(pout), np.array(mult)
##### original code ####
# p = asarray(p) * 1.0
# tol = abs(tol)
# p, indx = cmplx_sort(p)
# pout = []
# mult = []
# indx = -1
# curp = p[0] + 5 * tol
# sameroots = []
# for k in range(len(p)):
# tr = p[k]
# if abs(tr - curp) < tol:
# sameroots.append(tr)
# curp = comproot(sameroots)
# pout[indx] = curp
# mult[indx] += 1
# else:
# pout.append(tr)
# curp = tr
# sameroots = [tr]
# indx += 1
# mult.append(1)
# return array(pout), array(mult)
#------------------------------------------------------------------------------
def zplane(b=None, a=1, z=None, p=None, k=1, pn_eps=1e-3, analog=False,
plt_ax = None, plt_poles=True, style='square', anaCircleRad=0, lw=2,
mps = 10, mzs = 10, mpc = 'r', mzc = 'b', plabel = '', zlabel = ''):
"""
Plot the poles and zeros in the complex z-plane either from the
coefficients (`b,`a) of a discrete transfer function `H`(`z`) (zpk = False)
or directly from the zeros and poles (z,p) (zpk = True).
When only b is given, an FIR filter with all poles at the origin is assumed.
Parameters
----------
b : array_like
Numerator coefficients (transversal part of filter)
When b is not None, poles and zeros are determined from the coefficients
b and a
a : array_like (optional, default = 1 for FIR-filter)
Denominator coefficients (recursive part of filter)
z : array_like, default = None
Zeros
When b is None, poles and zeros are taken directly from z and p
p : array_like, default = None
Poles
analog : boolean (default: False)
When True, create a P/Z plot suitable for the s-plane, i.e. suppress
the unit circle (unless anaCircleRad > 0) and scale the plot for
a good display of all poles and zeros.
pn_eps : float (default : 1e-2)
Tolerance for separating close poles or zeros
plt_ax : handle to axes for plotting (default: None)
When no axes is specified, the current axes is determined via plt.gca()
plt_poles : Boolean (default : True)
Plot poles. This can be used to suppress poles for FIR systems
where all poles are at the origin.
style : string (default: 'square')
Style of the plot, for style == 'square' make scale of x- and y-
axis equal.
mps : integer (default: 10)
Size for pole marker
mzs : integer (default: 10)
Size for zero marker
mpc : char (default: 'r')
Pole marker colour
mzc : char (default: 'b')
Zero marker colour
lw : integer (default: 2)
Linewidth for unit circle
plabel, zlabel : string (default: '')
This string is passed to the plot command for poles and zeros and
can be displayed by legend()
Returns
-------
z, p, k : ndarray
Notes
-----
"""
# TODO:
# - polar option
# - add keywords for color of circle -> **kwargs
# - add option for multi-dimensional arrays and zpk data
# make sure that all inputs are arrays
b = np.atleast_1d(b)
a = np.atleast_1d(a)
z = np.atleast_1d(z) # make sure that p, z are arrays
p = np.atleast_1d(p)
if b.any(): # coefficients were specified
if len(b) < 2 and len(a) < 2:
logger.error('No proper filter coefficients: both b and a are scalars!')
return z, p, k
# The coefficients are less than 1, normalize the coefficients
if np.max(b) > 1:
kn = np.max(b)
b = b / float(kn)
else:
kn = 1.
if np.max(a) > 1:
kd = np.max(a)
a = a / abs(kd)
else:
kd = 1.
# Calculate the poles, zeros and scaling factor
p = np.roots(a)
z = np.roots(b)
k = kn/kd
elif not (len(p) or len(z)): # P/Z were specified
logger.error('Either b,a or z,p must be specified!')
return z, p, k
# find multiple poles and zeros and their multiplicities
if len(p) < 2: # single pole, [None] or [0]
if not p or p == 0: # only zeros, create equal number of poles at origin
p = np.array(0,ndmin=1) #
num_p = np.atleast_1d(len(z))
else:
num_p = [1.] # single pole != 0
else:
#p, num_p = sig.signaltools.unique_roots(p, tol = pn_eps, rtype='avg')
p, num_p = unique_roots(p, tol = pn_eps, rtype='avg')
# p = np.array(p); num_p = np.ones(len(p))
if len(z) > 0:
z, num_z = unique_roots(z, tol = pn_eps, rtype='avg')
# z = np.array(z); num_z = np.ones(len(z))
#z, num_z = sig.signaltools.unique_roots(z, tol = pn_eps, rtype='avg')
else:
num_z = []
ax = plt_ax#.subplot(111)
if analog == False:
# create the unit circle for the z-plane
uc = patches.Circle((0,0), radius=1, fill=False,
color='grey', ls='solid', zorder=1)
ax.add_patch(uc)
if style == 'square':
#r = 1.1
#ax.axis([-r, r, -r, r]) # overridden by next option
ax.axis('equal')
# ax.spines['left'].set_position('center')
# ax.spines['bottom'].set_position('center')
# ax.spines['right'].set_visible(True)
# ax.spines['top'].set_visible(True)
else: # s-plane
if anaCircleRad > 0:
# plot a circle with radius = anaCircleRad
uc = patches.Circle((0,0), radius=anaCircleRad, fill=False,
color='grey', ls='solid', zorder=1)
ax.add_patch(uc)
# plot real and imaginary axis
ax.axhline(lw=2, color = 'k', zorder=1)
ax.axvline(lw=2, color = 'k', zorder=1)
# Plot the zeros
ax.scatter(z.real, z.imag, s=mzs*mzs, zorder=2, marker = 'o',
facecolor = 'none', edgecolor = mzc, lw = lw, label=zlabel)
# and print their multiplicity
for i in range(len(z)):
if num_z[i] > 1:
ax.text(np.real(z[i]), np.imag(z[i]),' (' + str(num_z[i]) +')',
va = 'top', color=mzc)
if plt_poles:
# Plot the poles
ax.scatter(p.real, p.imag, s=mps*mps, zorder=2, marker='x',
color=mpc, lw=lw, label=plabel)
# and print their multiplicity
for i in range(len(p)):
if num_p[i] > 1:
ax.text(np.real(p[i]), np.imag(p[i]), ' (' + str(num_p[i]) +')',
va = 'bottom', color=mpc)
# =============================================================================
# # increase distance between ticks and labels
# # to give some room for poles and zeros
# for tick in ax.get_xaxis().get_major_ticks():
# tick.set_pad(12.)
# tick.label1 = tick._get_text1()
# for tick in ax.get_yaxis().get_major_ticks():
# tick.set_pad(12.)
# tick.label1 = tick._get_text1()
#
# =============================================================================
xl = ax.get_xlim(); Dx = max(abs(xl[1]-xl[0]), 0.05)
yl = ax.get_ylim(); Dy = max(abs(yl[1]-yl[0]), 0.05)
ax.set_xlim((xl[0]-Dx*0.05, max(xl[1]+Dx*0.05,0)))
ax.set_ylim((yl[0]-Dy*0.05, yl[1] + Dy*0.05))
return z, p, k
#------------------------------------------------------------------------------
def zplane_bak(b=None, a=1, z=None, p=None, k=1, pn_eps=1e-3, analog=False,
plt_ax = None, plt_poles=True, style='square', anaCircleRad=0, lw=2,
mps = 10, mzs = 10, mpc = 'r', mzc = 'b', plabel = '', zlabel = ''):
"""
Plot the poles and zeros in the complex z-plane either from the
coefficients (`b,`a) of a discrete transfer function `H`(`z`) (b specified)
or directly from the zeros and poles (z,p specified).
When only b is given, an FIR filter with all poles at the origin is assumed.
Parameters
----------
b : array_like
Numerator coefficients (transversal part of filter)
When b is not None, poles and zeros are determined from the coefficients
b and a
a : array_like (optional, default = 1 for FIR-filter)
Denominator coefficients (recursive part of filter)
z : array_like, default = None
Zeros
When b is None, poles and zeros are taken directly from z and p
p : array_like, default = None
Poles
analog : boolean (default: False)
When True, create a P/Z plot suitable for the s-plane, i.e. suppress
the unit circle (unless anaCircleRad > 0) and scale the plot for
a good display of all poles and zeros.
pn_eps : float (default : 1e-2)
Tolerance for separating close poles or zeros
plt_ax : handle to axes for plotting (default: None)
When no axes is specified, the current axes is determined via plt.gca()
plt_poles : Boolean (default : True)
Plot poles. This can be used to suppress poles for FIR systems
where all poles are at the origin.
style : string (default: 'square')
Style of the plot, for style == 'square' make scale of x- and y-
axis equal.
mps : integer (default: 10)
Size for pole marker
mzs : integer (default: 10)
Size for zero marker
mpc : char (default: 'r')
Pole marker colour
mzc : char (default: 'b')
Zero marker colour
lw : integer (default: 2)
Linewidth for unit circle
plabel, zlabel : string (default: '')
This string is passed to the plot command for poles and zeros and
can be displayed by legend()
Returns
-------
z, p, k : ndarray
Notes
-----
"""
# TODO:
# - polar option
# - add keywords for color of circle -> **kwargs
# - add option for multi-dimensional arrays and zpk data
# make sure that all inputs are arrays
b = np.atleast_1d(b)
a = np.atleast_1d(a)
z = np.atleast_1d(z) # make sure that p, z are arrays
p = np.atleast_1d(p)
if b.any(): # coefficients were specified
if len(b) < 2 and len(a) < 2:
logger.error('No proper filter coefficients: both b and a are scalars!')
return z, p, k
# The coefficients are less than 1, normalize the coefficients
if np.max(b) > 1:
kn = np.max(b)
b = b / float(kn)
else:
kn = 1.
if np.max(a) > 1:
kd = np.max(a)
a = a / abs(kd)
else:
kd = 1.
# Calculate the poles, zeros and scaling factor
p = np.roots(a)
z = np.roots(b)
k = kn/kd
elif not (len(p) or len(z)): # P/Z were specified
print('Either b,a or z,p must be specified!')
return z, p, k
# find multiple poles and zeros and their multiplicities
if len(p) < 2: # single pole, [None] or [0]
if not p or p == 0: # only zeros, create equal number of poles at origin
p = np.array(0,ndmin=1) #
num_p = np.atleast_1d(len(z))
else:
num_p = [1.] # single pole != 0
else:
#p, num_p = sig.signaltools.unique_roots(p, tol = pn_eps, rtype='avg')
p, num_p = unique_roots(p, tol = pn_eps, rtype='avg')
# p = np.array(p); num_p = np.ones(len(p))
if len(z) > 0:
z, num_z = unique_roots(z, tol = pn_eps, rtype='avg')
# z = np.array(z); num_z = np.ones(len(z))
#z, num_z = sig.signaltools.unique_roots(z, tol = pn_eps, rtype='avg')
else:
num_z = []
if not plt_ax:
ax = plt.gca()# fig.add_subplot(111)
else:
ax = plt_ax
if analog == False:
# create the unit circle for the z-plane
uc = patches.Circle((0,0), radius=1, fill=False,
color='grey', ls='solid', zorder=1)
ax.add_patch(uc)
if style == 'square':
r = 1.1
ax.axis([-r, r, -r, r])
ax.axis('equal')
# ax.spines['left'].set_position('center')
# ax.spines['bottom'].set_position('center')
# ax.spines['right'].set_visible(True)
# ax.spines['top'].set_visible(True)
else: # s-plane
if anaCircleRad > 0:
# plot a circle with radius = anaCircleRad
uc = patches.Circle((0,0), radius=anaCircleRad, fill=False,
color='grey', ls='solid', zorder=1)
ax.add_patch(uc)
# plot real and imaginary axis
ax.axhline(lw=2, color = 'k', zorder=1)
ax.axvline(lw=2, color = 'k', zorder=1)
# Plot the zeros
ax.scatter(z.real, z.imag, s=mzs*mzs, zorder=2, marker = 'o',
facecolor = 'none', edgecolor = mzc, lw = lw, label=zlabel)
# and print their multiplicity
for i in range(len(z)):
if num_z[i] > 1:
ax.text(np.real(z[i]), np.imag(z[i]),' (' + str(num_z[i]) +')',
va = 'top', color=mzc)
if plt_poles:
# Plot the poles
ax.scatter(p.real, p.imag, s=mps*mps, zorder=2, marker='x',
color=mpc, lw=lw, label=plabel)
# and print their multiplicity
for i in range(len(p)):
if num_p[i] > 1:
ax.text(np.real(p[i]), np.imag(p[i]), ' (' + str(num_p[i]) +')',
va = 'bottom', color=mpc)
# increase distance between ticks and labels
# to give some room for poles and zeros
for tick in ax.get_xaxis().get_major_ticks():
tick.set_pad(12.)
tick.label1 = tick._get_text1()
for tick in ax.get_yaxis().get_major_ticks():
tick.set_pad(12.)
tick.label1 = tick._get_text1()
xl = ax.get_xlim(); Dx = max(abs(xl[1]-xl[0]), 0.05)
yl = ax.get_ylim(); Dy = max(abs(yl[1]-yl[0]), 0.05)
ax.set_xlim((xl[0]-Dx*0.05, max(xl[1]+Dx*0.05,0)))
ax.set_ylim((yl[0]-Dy*0.05, yl[1] + Dy*0.05))
return z, p, k
#------------------------------------------------------------------------------
def impz(b, a=1, FS=1, N=1, step=False):
"""
Calculate impulse response of a discrete time filter, specified by
numerator coefficients b and denominator coefficients a of the system
function H(z).
When only b is given, the impulse response of the transversal (FIR)
filter specified by b is calculated.
Parameters
----------
b : array_like
Numerator coefficients (transversal part of filter)
a : array_like (optional, default = 1 for FIR-filter)
Denominator coefficients (recursive part of filter)
FS : float (optional, default: FS = 1)
Sampling frequency.
N : float (optional, default N=1 for automatic calculation)
Number of calculated points.
Default: N = len(b) for FIR filters, N = 100 for IIR filters
step: boolean (optional, default: step=False)
plot step response instead of impulse response
Returns
-------
hn : ndarray with length N (see above)
td : ndarray containing the time steps with same
Examples
--------
>>> b = [1,2,3] # Coefficients of H(z) = 1 + 2 z^2 + 3 z^3
>>> h, n = dsp_lib.impz(b)
"""
try: len(a) #len_a = len(a)
except TypeError:
# a has len = 1 -> FIR-Filter
impulse = np.repeat(0.,len(b)) # create float array filled with 0.
try: len(b)
except TypeError:
print('No proper filter coefficients: len(a) = len(b) = 1 !')
else:
try: len(b)
except TypeError: b = [b,] # convert scalar to array with len = 1
impulse = np.repeat(0.,100) # IIR-Filter
if N > 1:
impulse = np.repeat(0.,N)
impulse[0] =1.0 # create dirac impulse
hn = np.array(sig.lfilter(b,a,impulse)) # calculate impulse response
td = np.arange(len(hn)) / FS
if step:
hn = np.cumsum(hn) # integrate impulse response to get step response
return hn, td
#==================================================================
def grpdelay(b, a=1, nfft=512, whole='none', analog=False, Fs=2.*pi):
#==================================================================
"""
Calculate group delay of a discrete time filter, specified by
numerator coefficients `b` and denominator coefficients `a` of the system
function `H` ( `z`).
When only `b` is given, the group delay of the transversal (FIR)
filter specified by `b` is calculated.
Parameters
----------
b : array_like
Numerator coefficients (transversal part of filter)
a : array_like (optional, default = 1 for FIR-filter)
Denominator coefficients (recursive part of filter)
whole : string (optional, default : 'none')
Only when whole = 'whole' calculate group delay around
the complete unit circle (0 ... 2 pi)
N : integer (optional, default: 512)
Number of FFT-points
FS : float (optional, default: FS = 2*pi)
Sampling frequency.
Returns
-------
tau_g : ndarray
The group delay
w : ndarray
The angular frequency points where the group delay was computed
Notes
-----
The group delay :math:`\\tau_g(\\omega)` of discrete and continuous time
systems is defined by
.. math::
\\tau_g(\\omega) = - \\phi'(\\omega)
= -\\frac{\\partial \\phi(\\omega)}{\\partial \\omega}
= -\\frac{\\partial }{\\partial \\omega}\\angle H( \\omega)
A useful form for calculating the group delay is obtained by deriving the
*logarithmic* frequency response in polar form as described in [JOS]_ for
discrete time systems:
.. math::
\\ln ( H( \\omega))
= \\ln \\left({H_A( \\omega)} e^{j \\phi(\\omega)} \\right)
= \\ln \\left({H_A( \\omega)} \\right) + j \\phi(\\omega)
\\Rightarrow \\; \\frac{\\partial }{\\partial \\omega} \\ln ( H( \\omega))
= \\frac{H_A'( \\omega)}{H_A( \\omega)} + j \\phi'(\\omega)
where :math:`H_A(\\omega)` is the amplitude response. :math:`H_A(\\omega)` and
its derivative :math:`H_A'(\\omega)` are real-valued, therefore, the group
delay can be calculated from
.. math::
\\tau_g(\\omega) = -\\phi'(\\omega) =
-\\Im \\left\\{ \\frac{\\partial }{\\partial \\omega}
\\ln ( H( \\omega)) \\right\\}
=-\\Im \\left\\{ \\frac{H'(\\omega)}{H(\\omega)} \\right\\}
The derivative of a polynome :math:`P(s)` (continuous-time system) or :math:`P(z)`
(discrete-time system) w.r.t. :math:`\\omega` is calculated by:
.. math::
\\frac{\\partial }{\\partial \\omega} P(s = j \\omega)
= \\frac{\\partial }{\\partial \\omega} \\sum_{k = 0}^N c_k (j \\omega)^k
= j \\sum_{k = 0}^{N-1} (k+1) c_{k+1} (j \\omega)^{k}
= j P_R(s = j \\omega)
\\frac{\\partial }{\\partial \\omega} P(z = e^{j \\omega T})
= \\frac{\\partial }{\\partial \\omega} \\sum_{k = 0}^N c_k e^{-j k \\omega T}
= -jT \\sum_{k = 0}^{N} k c_{k} e^{-j k \\omega T}
= -jT P_R(z = e^{j \\omega T})
where :math:`P_R` is the "ramped" polynome, i.e. its `k` th coefficient is
multiplied by `k` resp. `k` + 1.
yielding:
.. math::
\\tau_g(\\omega) = -\\Im \\left\\{ \\frac{H'(\\omega)}{H(\\omega)} \\right\\}
\\quad \\text{ resp. } \\quad
\\tau_g(\\omega) = -\\Im \\left\\{ \\frac{H'(e^{j \\omega T})}
{H(e^{j \\omega T})} \\right\\}
where::
(H'(e^jwT)) ( H_R(e^jwT)) (H_R(e^jwT))
tau_g(w) = -im |---------| = -im |-jT ----------| = T re |----------|
( H(e^jwT)) ( H(e^jwT) ) ( H(e^jwT) )
where :math:`H(e^{j\\omega T})` is calculated via the DFT at NFFT points and
the derivative
of the polynomial terms :math:`b_k z^-k` using :math:`\\partial / \\partial w b_k e^-jkwT` = -b_k jkT e^-jkwT.
This is equivalent to muliplying the polynome with a ramp `k`,
yielding the "ramped" function H_R(e^jwT).
For analog functions with b_k s^k the procedure is analogous, but there is no
sampling time and the exponent is positive.
.. [JOS] Julius O. Smith III, "Numerical Computation of Group Delay" in
"Introduction to Digital Filters with Audio Applications",
Center for Computer Research in Music and Acoustics (CCRMA),
Stanford University, http://ccrma.stanford.edu/~jos/filters/Numerical_Computation_Group_Delay.html, referenced 2014-04-02,
.. [Lyons] Richard Lyons, "Understanding Digital Signal Processing", 3rd Ed.,
Prentice Hall, 2010.
Examples
--------
>>> b = [1,2,3] # Coefficients of H(z) = 1 + 2 z^2 + 3 z^3
>>> tau_g, td = dsp_lib.grpdelay(b)
"""
## If the denominator of the computation becomes too small, the group delay
## is set to zero. (The group delay approaches infinity when
## there are poles or zeros very close to the unit circle in the z plane.)
##
## Theory: group delay, g(w) = -d/dw [arg{H(e^jw)}], is the rate of change of
## phase with respect to frequency. It can be computed as:
##
## d/dw H(e^-jw)
## g(w) = -------------
## H(e^-jw)
##
## where
## H(z) = B(z)/A(z) = sum(b_k z^k)/sum(a_k z^k).
##
## By the quotient rule,
## A(z) d/dw B(z) - B(z) d/dw A(z)
## d/dw H(z) = -------------------------------
## A(z) A(z)
## Substituting into the expression above yields:
## A dB - B dA
## g(w) = ----------- = dB/B - dA/A
## A B
##
## Note that,
## d/dw B(e^-jw) = sum(k b_k e^-jwk)
## d/dw A(e^-jw) = sum(k a_k e^-jwk)
## which is just the FFT of the coefficients multiplied by a ramp.
##
## As a further optimization when nfft>>length(a), the IIR filter (b,a)
## is converted to the FIR filter conv(b,fliplr(conj(a))).
if whole !='whole':
nfft = 2*nfft
nfft = int(nfft)
#
w = Fs * np.arange(0, nfft)/nfft # create frequency vector
try: len(a)
except TypeError:
a = 1; oa = 0 # a is a scalar or empty -> order of a = 0
c = b
try: len(b)
except TypeError: print('No proper filter coefficients: len(a) = len(b) = 1 !')
else:
oa = len(a)-1 # order of denom. a(z) resp. a(s)
c = np.convolve(b,a[::-1]) # a[::-1] reverses denominator coeffs a
# c(z) = b(z) * a(1/z)*z^(-oa)
try: len(b)
except TypeError: b=1; ob=0 # b is a scalar or empty -> order of b = 0
else:
ob = len(b)-1 # order of b(z)
if analog:
a_b = np.convolve(a,b)
if ob > 1:
br_a = np.convolve(b[1:] * np.arange(1,ob), a)
else:
br_a = 0
ar_b = np.convolve(a[1:] * np.arange(1,oa), b)
num = np.fft.fft(ar_b - br_a, nfft)
den = np.fft.fft(a_b,nfft)
else:
oc = oa + ob # order of c(z)
cr = c * np.arange(0,oc+1) # multiply with ramp -> derivative of c wrt 1/z
num = np.fft.fft(cr,nfft) #
den = np.fft.fft(c,nfft) #
#
minmag = 10. * np.spacing(1) # equivalent to matlab "eps"
polebins = np.where(abs(den) < minmag)[0] # find zeros of denominator
# polebins = np.where(abs(num) < minmag)[0] # find zeros of numerator
if np.size(polebins) > 0: # check whether polebins array is empty
print('*** grpdelay warning: group delay singular -> setting to 0 at:')
for i in polebins:
print ('f = {0} '.format((Fs*i/nfft)))
num[i] = 0
den[i] = 1
if analog:
tau_g = np.real(num / den)
else:
tau_g = np.real(num / den) - oa
#
if whole !='whole':
nfft = nfft//2
tau_g = tau_g[0:nfft]
w = w[0:nfft]
return tau_g, w
def grp_delay_ana(b, a, w):
"""
Calculate the group delay of an anlog filter.
"""
w, H = sig.freqs(b, a, w)
H_angle = np.unwrap(np.angle(H))
# tau_g = np.zeros(len(w)-1)
tau_g = (H_angle[1:]-H_angle[:-1])/(w[0]-w[1])
return tau_g, w[:-1]
#==================================================================
def format_ticks(xy, scale, format="%.1f"):
#==================================================================
"""
Reformat numbers at x or y - axis. The scale can be changed to display
e.g. MHz instead of Hz. The number format can be changed as well.
Parameters
----------
xy : string, either 'x', 'y' or 'xy'
select corresponding axis (axes) for reformatting
scale : real,
format : string,
define C-style number formats
Returns
-------
nothing
Examples
--------
>>> format_ticks('x',1000.)
Scales all numbers of x-Axis by 1000, e.g. for displaying ms instead of s.
>>> format_ticks('xy',1., format = "%.2f")
Two decimal places for numbers on x- and y-axis
"""
if xy == 'x' or xy == 'xy':
locx,labelx = plt.xticks() # get location and content of xticks
plt.xticks(locx, map(lambda x: format % x, locx*scale))
if xy == 'y' or xy == 'xy':
locy,labely = plt.yticks() # get location and content of xticks
plt.yticks(locy, map(lambda y: format % y, locy*scale))
#==================================================================
def lim_eps(a, eps):
#==================================================================
"""
Return min / max of an array a, increased by eps*(max(a) - min(a)).
Handy for nice looking axes labeling.
"""
mylim = (min(a) - (max(a)-min(a))*eps, max(a) + (max(a)-min(a))*eps)
return mylim
#========================================================
abs = absolute
def oddround(x):
"""Return the nearest odd integer from x."""
return x-mod(x,2)+1
def oddceil(x):
"""Return the smallest odd integer not less than x."""
return oddround(x+1)
def remlplen_herrmann(fp,fs,dp,ds):
"""
Determine the length of the low pass filter with passband frequency
fp, stopband frequency fs, passband ripple dp, and stopband ripple ds.
fp and fs must be normalized with respect to the sampling frequency.
Note that the filter order is one less than the filter length.
Uses approximation algorithm described by Herrmann et al.:
O. Herrmann, L.R. Raviner, and D.S.K. Chan, Practical Design Rules for
Optimum Finite Impulse Response Low-Pass Digital Filters, Bell Syst. Tech.
Jour., 52(6):769-799, Jul./Aug. 1973.
"""
dF = fs-fp
a = [5.309e-3,7.114e-2,-4.761e-1,-2.66e-3,-5.941e-1,-4.278e-1]
b = [11.01217, 0.51244]
Dinf = log10(ds)*(a[0]*log10(dp)**2+a[1]*log10(dp)+a[2])+ \
a[3]*log10(dp)**2+a[4]*log10(dp)+a[5]
f = b[0]+b[1]*(log10(dp)-log10(ds))
N1 = Dinf/dF-f*dF+1
return int(oddround(N1))
def remlplen_kaiser(fp,fs,dp,ds):
"""
Determine the length of the low pass filter with passband frequency
fp, stopband frequency fs, passband ripple dp, and stopband ripple ds.
fp and fs must be normalized with respect to the sampling frequency.
Note that the filter order is one less than the filter length.
Uses approximation algorithm described by Kaiser:
J.F. Kaiser, Nonrecursive Digital Filter Design Using I_0-sinh Window
function, Proc. IEEE Int. Symp. Circuits and Systems, 20-23, April 1974.
"""
dF = fs-fp
N2 = (-20*log10(sqrt(dp*ds))-13.0)/(14.6*dF)+1.0
return int(oddceil(N2))
def remlplen_ichige(fp,fs,dp,ds):
"""
Determine the length of the low pass filter with passband frequency
fp, stopband frequency fs, passband ripple dp, and stopband ripple ds.
fp and fs must be normalized with respect to the sampling frequency.
Note that the filter order is one less than the filter length.
Uses approximation algorithm described by Ichige et al.:
K. Ichige, M. Iwaki, and R. Ishii, Accurate Estimation of Minimum
Filter Length for Optimum FIR Digital Filters, IEEE Transactions on
Circuits and Systems, 47(10):1008-1017, October 2000.
"""
dF = fs-fp
v = lambda dF,dp:2.325*((-log10(dp))**-0.445)*dF**(-1.39)
g = lambda fp,dF,d:(2.0/pi)*arctan(v(dF,dp)*(1.0/fp-1.0/(0.5-dF)))
h = lambda fp,dF,c:(2.0/pi)*arctan((c/dF)*(1.0/fp-1.0/(0.5-dF)))
Nc = ceil(1.0+(1.101/dF)*(-log10(2.0*dp))**1.1)
Nm = (0.52/dF)*log10(dp/ds)*(-log10(dp))**0.17
N3 = ceil(Nc*(g(fp,dF,dp)+g(0.5-dF-fp,dF,dp)+1.0)/3.0)
DN = ceil(Nm*(h(fp,dF,1.1)-(h(0.5-dF-fp,dF,0.29)-1.0)/2.0))
N4 = N3+DN
return int(N4)
def remezord(freqs,amps,rips,Hz=1,alg='ichige'):
"""
Filter parameter selection for the Remez exchange algorithm.
Calculate the parameters required by the Remez exchange algorithm to
construct a finite impulse response (FIR) filter that approximately
meets the specified design.
Parameters
----------
freqs : list
A monotonic sequence of band edges specified in Hertz. All elements
must be non-negative and less than 1/2 the sampling frequency as
given by the Hz parameter. The band edges "0" and "f_S / 2" do not
have to be specified, hence 2 * number(amps) - 2 freqs are needed.
amps : list
A sequence containing the amplitudes of the signal to be
filtered over the various bands, e.g. 1 for the passband, 0 for the
stopband and 0.42 for some intermediate band.
rips : list
A list with the peak ripples (linear, not in dB!) for each band. For
the stop band this is equivalent to the minimum attenuation.
Hz : float
Sampling frequency
alg : string
Filter length approximation algorithm. May be either 'herrmann',
'kaiser' or 'ichige'. Depending on the specifications, some of
the algorithms may give better results than the others.
Returns
-------
numtaps,bands,desired,weight -- See help for the remez function.
Examples
--------
We want to design a lowpass with the band edges of 40 resp. 50 Hz and a
sampling frequency of 200 Hz, a passband peak ripple of 10%
and a stop band ripple of 0.01 or 40 dB.
>>> (L, F, A, W) = dsp.remezord([40, 50], [1, 0], [0.1, 0.01], Hz = 200)
Notes
-----
Supplies remezord method according to Scipy Ticket #475
http://projects.scipy.org/scipy/ticket/475
https://github.com/scipy/scipy/issues/1002
https://github.com/thorstenkranz/eegpy/blob/master/eegpy/filter/remezord.py
"""
# Make sure the parameters are floating point numpy arrays:
freqs = asarray(freqs,'d')
amps = asarray(amps,'d')
rips = asarray(rips,'d')
# Scale ripples with respect to band amplitudes:
rips /= (amps+(amps==0.0))
# Normalize input frequencies with respect to sampling frequency:
freqs /= Hz
# Select filter length approximation algorithm:
if alg == 'herrmann':
remlplen = remlplen_herrmann
elif alg == 'kaiser':
remlplen = remlplen_kaiser
elif alg == 'ichige':
remlplen = remlplen_ichige
else:
raise ValueError('Unknown filter length approximation algorithm.')
# Validate inputs:
if any(freqs > 0.5):
raise ValueError('Frequency band edges must not exceed the Nyquist frequency.')
if any(freqs < 0.0):
raise ValueError('Frequency band edges must be nonnegative.')
if any(rips <= 0.0):
raise ValueError('Ripples must be nonnegative and non-zero.')
if len(amps) != len(rips):
raise ValueError('Number of amplitudes must equal number of ripples.')
if len(freqs) != 2*(len(amps)-1):
raise ValueError('Number of band edges must equal 2*(number of amplitudes-1)')
# Find the longest filter length needed to implement any of the
# low-pass or high-pass filters with the specified edges:
f1 = freqs[0:-1:2]
f2 = freqs[1::2]
L = 0
for i in range(len(amps)-1):
L = max((L,
remlplen(f1[i],f2[i],rips[i],rips[i+1]),
remlplen(0.5-f2[i],0.5-f1[i],rips[i+1],rips[i])))
# Cap the sequence of band edges with the limits of the digital frequency
# range:
bands = hstack((0.0,freqs,0.5))
# The filter design weights correspond to the ratios between the maximum
# ripple and all of the other ripples:
weight = max(rips)/rips
return [L,bands,amps,weight]
#######################################
# If called directly, do some example #
#######################################
if __name__=='__main__':
pass
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment