diff --git a/applications/lofar2/model/pfb_os/dsp_fpga_lib.py b/applications/lofar2/model/pfb_os/dsp_fpga_lib.py
index 55077a0538b8fa06f894d7b6861b04484d2c9019..629fd59d160f6be6cac9d7ce8a06f3bbd60266fa 100644
--- a/applications/lofar2/model/pfb_os/dsp_fpga_lib.py
+++ b/applications/lofar2/model/pfb_os/dsp_fpga_lib.py
@@ -1,4 +1,5 @@
# -*- coding: utf-8 -*-
+# flake8: noqa
"""
Created on Mon Apr 30 10:29:42 2012
dsp_fpga_lib (Version 8)
@@ -6,18 +7,18 @@ dsp_fpga_lib (Version 8)
from: https://github.com/chipmuenk/dsp
"""
import sys
-import string # needed for remezord?
+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
+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
+from matplotlib import patches
__version__ = "0.9"
def versions():
@@ -26,20 +27,21 @@ def versions():
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.family' : 'sans-serif', #'serif',
'font.style' : 'normal',
- #'mathtext.fontset' : 'stixsans',#'stix',
- #'mathtext.fallback_to_cm' : True,
+ # '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.linewidth' : 1, # linewidth for coordinate system
+ 'axes.formatter.use_mathtext': True, # use mathtext for scientific notation.
#
'axes.grid' : True,
'grid.color' : '#222222',
@@ -55,7 +57,7 @@ mpl_rc = {'lines.linewidth' : 1.5,
'xtick.color' : 'black',
'ytick.color' : 'black',
#
- 'figure.figsize' : (7,4), # default figure size in inches
+ 'figure.figsize' : (7,4), # default figure size in inches
'figure.facecolor' : 'white',
'figure.edgecolor' : '#808080',
'figure.dpi' : 100,
@@ -64,20 +66,21 @@ mpl_rc = {'lines.linewidth' : 1.5,
'savefig.edgecolor' : 'white',
'savefig.bbox' : 'tight',
'savefig.pad_inches' : 0,
- 'hatch.color' : '#808080', # mpl >= 2.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
+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
+ 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
+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
@@ -86,19 +89,19 @@ def H_mag(zaehler, nenner, z, lim):
# limvec = lim * np.ones(len(z))
try: len(zaehler)
except TypeError:
- z_val = abs(zaehler) # zaehler is a scalar
+ z_val = abs(zaehler) # zaehler is a scalar
else:
- z_val = abs(np.polyval(zaehler,z)) # evaluate zaehler at z
+ z_val = abs(np.polyval(zaehler, z)) # evaluate zaehler at z
try: len(nenner)
except TypeError:
- n_val = nenner # nenner is a scalar
+ n_val = nenner # nenner is a scalar
else:
- n_val = abs(np.polyval(nenner,z))
+ n_val = abs(np.polyval(nenner, z))
- return np.minimum((z_val/n_val),lim)
+ return np.minimum((z_val/n_val), lim)
-#----------------------------------------------
+# ----------------------------------------------
# from scipy.sig.signaltools.py:
def cmplx_sort(p):
"sort roots based on magnitude."
@@ -109,76 +112,73 @@ def cmplx_sort(p):
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')
+ """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.
- def manhattan(a,b):
- """
- Manhattan distance between a and b
- """
+ 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
- """
+
+ 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']:
@@ -196,41 +196,37 @@ uniq, mult = sp.signal.unique_roots(vals, rtype='avg')
else:
raise TypeError(rdist)
- mult = [] # initialize list for multiplicities
- pout = [] # initialize list for reduced output list of roots
+ 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
+ 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?
+ 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
+ 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"
+ 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
-
+ pout.append(comproot(sameroots)) # avg/mean/max of mult. root
else:
- p,indx = cmplx_sort(p)
+ p, indx = cmplx_sort(p)
indx = -1
- curp = p[0] + 5 * tol # needed to avoid "self-detection" ?
+ 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:
@@ -246,7 +242,6 @@ uniq, mult = sp.signal.unique_roots(vals, rtype='avg')
sameroots = [tr]
indx += 1
mult.append(1)
-
return np.array(pout), np.array(mult)
##### original code ####
@@ -275,10 +270,9 @@ uniq, mult = sp.signal.unique_roots(vals, rtype='avg')
#------------------------------------------------------------------------------
-
-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 = ''):
+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)
@@ -341,12 +335,10 @@ def zplane(b=None, a=1, z=None, p=None, k=1, pn_eps=1e-3, analog=False,
This string is passed to the plot command for poles and zeros and
can be displayed by legend()
-
Returns
-------
z, p, k : ndarray
-
Notes
-----
"""
@@ -358,10 +350,10 @@ def zplane(b=None, a=1, z=None, p=None, k=1, pn_eps=1e-3, analog=False,
# 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
+ z = np.atleast_1d(z) # make sure that p, z are arrays
p = np.atleast_1d(p)
- if b.any(): # coefficients were specified
+ 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
@@ -383,97 +375,97 @@ def zplane(b=None, a=1, z=None, p=None, k=1, pn_eps=1e-3, analog=False,
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!')
+ 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) #
+ 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
+ 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, 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, 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')
+ # 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,
+ if plt_ax:
+ 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)
- # 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))
+ 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).
+ 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.
@@ -532,12 +524,10 @@ def zplane_bak(b=None, a=1, z=None, p=None, k=1, pn_eps=1e-3, analog=False,
This string is passed to the plot command for poles and zeros and
can be displayed by legend()
-
Returns
-------
z, p, k : ndarray
-
Notes
-----
"""
@@ -549,10 +539,10 @@ def zplane_bak(b=None, a=1, z=None, p=None, k=1, pn_eps=1e-3, analog=False,
# 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
+ z = np.atleast_1d(z) # make sure that p, z are arrays
p = np.atleast_1d(p)
- if b.any(): # coefficients were specified
+ 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
@@ -574,36 +564,36 @@ def zplane_bak(b=None, a=1, z=None, p=None, k=1, pn_eps=1e-3, analog=False,
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!')
+ 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
+ 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
+ num_p = [1.] # single pole != 0
else:
- #p, num_p = sig.signaltools.unique_roots(p, tol = pn_eps, rtype='avg')
+ # 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')
+ # 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)
+ 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,
+ uc = patches.Circle((0, 0), radius=1, fill=False,
color='grey', ls='solid', zorder=1)
ax.add_patch(uc)
if style == 'square':
@@ -614,11 +604,10 @@ def zplane_bak(b=None, a=1, z=None, p=None, k=1, pn_eps=1e-3, analog=False,
# 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,
+ uc = patches.Circle((0, 0), radius=anaCircleRad, fill=False,
color='grey', ls='solid', zorder=1)
ax.add_patch(uc)
# plot real and imaginary axis
@@ -626,13 +615,13 @@ def zplane_bak(b=None, a=1, z=None, p=None, k=1, pn_eps=1e-3, analog=False,
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)
+ 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)
+ 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',
@@ -641,7 +630,7 @@ def zplane_bak(b=None, a=1, z=None, p=None, k=1, pn_eps=1e-3, analog=False,
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)
+ va='bottom', color=mpc)
# increase distance between ticks and labels
# to give some room for poles and zeros
@@ -652,13 +641,15 @@ def zplane_bak(b=None, a=1, z=None, p=None, k=1, pn_eps=1e-3, analog=False,
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))
-
+ 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):
"""
@@ -692,33 +683,33 @@ def impz(b, a=1, FS=1, N=1, step=False):
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
+ >>> 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.
+ 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
+ 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
+ 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
+ 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):
#==================================================================
@@ -748,13 +739,11 @@ def grpdelay(b, a=1, nfft=512, whole='none', analog=False, Fs=2.*pi):
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
@@ -820,7 +809,6 @@ def grpdelay(b, a=1, nfft=512, whole='none', analog=False, Fs=2.*pi):
\\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))
@@ -833,13 +821,9 @@ def grpdelay(b, a=1, nfft=512, whole='none', analog=False, Fs=2.*pi):
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),
@@ -850,82 +834,81 @@ def grpdelay(b, a=1, nfft=512, whole='none', analog=False, Fs=2.*pi):
Examples
--------
- >>> b = [1,2,3] # Coefficients of H(z) = 1 + 2 z^2 + 3 z^3
+ >>> 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 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 = 2 * nfft
nfft = int(nfft)
-#
- w = Fs * np.arange(0, nfft)/nfft # create frequency vector
+
+ 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
+ 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 = 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
+ 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)
+ ob = len(b) - 1 # order of b(z)
if analog:
- a_b = np.convolve(a,b)
+ a_b = np.convolve(a, b)
if ob > 1:
- br_a = np.convolve(b[1:] * np.arange(1,ob), a)
+ 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)
+ 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)
+ 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
+ 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
+ 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:
@@ -937,24 +920,24 @@ def grpdelay(b, a=1, nfft=512, whole='none', analog=False, Fs=2.*pi):
tau_g = np.real(num / den)
else:
tau_g = np.real(num / den) - oa
-#
- if whole !='whole':
- nfft = nfft//2
+
+ 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.
- """
+ """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])
+ tau_g = (H_angle[1:] - H_angle[:-1]) / (w[0] - w[1])
return tau_g, w[:-1]
+
#==================================================================
def format_ticks(xy, scale, format="%.1f"):
#==================================================================
@@ -976,20 +959,20 @@ def format_ticks(xy, scale, format="%.1f"):
-------
nothing
-
Examples
--------
- >>> format_ticks('x',1000.)
+ >>> 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")
+ >>> 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))
+ 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))
+ 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):
@@ -998,26 +981,26 @@ 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)
+ 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
- return x-mod(x,2)+1
def oddceil(x):
"""Return the smallest odd integer not less than x."""
+ return oddround(x + 1)
- return oddround(x+1)
-def remlplen_herrmann(fp,fs,dp,ds):
+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.
@@ -1030,18 +1013,17 @@ def remlplen_herrmann(fp,fs,dp,ds):
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]
+ 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
-
+ 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):
+
+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.
@@ -1053,13 +1035,12 @@ def remlplen_kaiser(fp,fs,dp,ds):
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
-
+ 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):
+
+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.
@@ -1070,22 +1051,20 @@ def remlplen_ichige(fp,fs,dp,ds):
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
-
+ 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.
+
+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
@@ -1131,20 +1110,18 @@ def remezord(freqs,amps,rips,Hz=1,alg='ichige'):
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')
+ freqs = asarray(freqs, 'd')
+ amps = asarray(amps, 'd')
+ rips = asarray(rips, 'd')
# Scale ripples with respect to band amplitudes:
- rips /= (amps+(amps==0.0))
+ rips /= (amps + (amps==0.0))
# Normalize input frequencies with respect to sampling frequency:
freqs /= Hz
@@ -1168,10 +1145,9 @@ def remezord(freqs,amps,rips,Hz=1,alg='ichige'):
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):
+ 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]
@@ -1179,18 +1155,18 @@ def remezord(freqs,amps,rips,Hz=1,alg='ichige'):
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])))
+ 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))
+ 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
+ weight = max(rips) / rips
+ return [L, bands, amps, weight]
- return [L,bands,amps,weight]
#######################################
# If called directly, do some example #