Python routine for conversion of IIR filter into cascade form
wzab
# This is a shell archive (produced by GNU sharutils 4.6.3).
# To extract the files from this archive, save it to some FILE, remove
# everything before the `#!/bin/sh' line above, then type `sh FILE'.
#
lock_dir=_sh08680
# Made on 2009-05-09 15:30 CEST by <wzab@wzab>.
# Source directory was `/tmp/iir'.
#
# Existing files will *not* be overwritten, unless `-c' is specified.
#
# This shar contains:
# length mode name
# ------ ---------- ------------------------------------------
# 8460 -rw------- iir2casc.py
#
MD5SUM=${MD5SUM-md5sum}
f=`${MD5SUM} --version | egrep '^md5sum .*(core|text)utils'`
test -n "${f}" && md5check=true || md5check=false
${md5check} || \
echo 'Note: not verifying md5sums. Consider installing GNU coreutils.'
save_IFS="${IFS}"
IFS="${IFS}:"
gettext_dir=FAILED
locale_dir=FAILED
first_param="$1"
for dir in $PATH
do
if test "$gettext_dir" = FAILED && test -f $dir/gettext \
&& ($dir/gettext --version >/dev/null 2>&1)
then
case `$dir/gettext --version 2>&1 | sed 1q` in
*GNU*) gettext_dir=$dir ;;
esac
fi
if test "$locale_dir" = FAILED && test -f $dir/shar \
&& ($dir/shar --print-text-domain-dir >/dev/null 2>&1)
then
locale_dir=`$dir/shar --print-text-domain-dir`
fi
done
IFS="$save_IFS"
if test "$locale_dir" = FAILED || test "$gettext_dir" = FAILED
then
echo=echo
else
TEXTDOMAINDIR=$locale_dir
export TEXTDOMAINDIR
TEXTDOMAIN=sharutils
export TEXTDOMAIN
echo="$gettext_dir/gettext -s"
fi
if (echo "testing\c"; echo 1,2,3) | grep c >/dev/null
then if (echo -n test; echo 1,2,3) | grep n >/dev/null
then shar_n= shar_c='
'
else shar_n=-n shar_c= ; fi
else shar_n= shar_c='\c' ; fi
f=shar-touch.$$
st1=200112312359.59
st2=123123592001.59
st2tr=123123592001.5 # old SysV 14-char limit
st3=1231235901
if touch -am -t ${st1} ${f} >/dev/null 2>&1 && \
test ! -f ${st1} && test -f ${f}; then
shar_touch='touch -am -t $1$2$3$4$5$6.$7 "$8"'
elif touch -am ${st2} ${f} >/dev/null 2>&1 && \
test ! -f ${st2} && test ! -f ${st2tr} && test -f ${f}; then
shar_touch='touch -am $3$4$5$6$1$2.$7 "$8"'
elif touch -am ${st3} ${f} >/dev/null 2>&1 && \
test ! -f ${st3} && test -f ${f}; then
shar_touch='touch -am $3$4$5$6$2 "$8"'
else
shar_touch=:
echo
${echo} 'WARNING: not restoring timestamps. Consider getting and'
${echo} 'installing GNU `touch'\'', distributed in GNU coreutils...'
echo
fi
rm -f ${st1} ${st2} ${st2tr} ${st3} ${f}
#
if test ! -d ${lock_dir}
then : ; else ${echo} 'lock directory '${lock_dir}' exists'
exit 1
fi
if mkdir ${lock_dir}
then ${echo} 'x - created lock directory `'${lock_dir}\''.'
else ${echo} 'x - failed to create lock directory `'${lock_dir}\''.'
exit 1
fi
# ============= iir2casc.py ==============
if test -f 'iir2casc.py' && test "$first_param" != -c; then
${echo} 'x -SKIPPING iir2casc.py (file already exists)'
else
${echo} 'x - extracting iir2casc.py (text)'
sed 's/^X//' << 'SHAR_EOF' > 'iir2casc.py' &&
#!/usr/bin/python
X
"""
This is the first attempt to prepare a function for conversion
of IIR filter in the direct form into a cascade form.
It requires a lot of testing and polishing yet...
Any suggestions are appreciated.
X
The data representation is probably not optimal.
X
I heavily rely on list append and pop operations,
so the code is probably not very fast :-(
X
This is PUBLIC DOMAIN code, written by Wojciech M. Zabolotny
wzab<at>ise.pw.edu.pl 6-9.05.2009
"""
import scipy
import scipy.signal
X
class sos:
X """
X Second order section (in some cases it may be also
X the first order section!!!)
X Fields:
X - num - numerator
X - den - denominator
X """
X def __init__(self,n,d):
X self.num=n
X self.den=d
X def __repr__(self):
X l="n:"+str(self.num.c)+" d:"+str(self.den.c)
X return l
X
X
class sos_n_or_d:
X """
X Second order numerator or denominator
X """
X def __init__(self,r,c):
X self.r=r
X self.c=c
X
class paired_roots:
X """
X Objects storing mainly the paired roots
X - complex (in "conjs" list)
X - real opposite roots (in "real_pairs" list)
X Additionally single real roots are kept in the "reals" list
X """
X def __init__(self):
X """
X Initialize the object, create empty lists: reals, conjs and real_pairs
X """
X self.reals=[]
X self.real_pairs=[]
X self.conjs=[]
X def add_real(self,r):
X """
X Add single real root
X """
X self.reals.append(r)
X def add_real_pair(self,r1,r2):
X """
X Add the pair of opposite real roots
X The positive root should be the first in the tupple
X """
X if r1>r2:
X self.real_pairs.append((r1,r2))
X else:
X self.real_pairs.append((r2,r1))
X def add_conjs(self,c1,c2):
X """
X Add the pair of conjugate roots.
X The root with imaginary part > 0 should be the first in the tupple
X """
X if scipy.imag(c1)>scipy.imag(c2):
X self.conjs.append((c1,c2))
X else:
X self.conjs.append((c2,c1))
X def get_sos_n_or_d(self):
X """
X This function returns the pair of roots
X and the coefficients of the second order polynomial,
X which may be used as the numerator or denominator
X of the second order section (class sos_n_or_d)
X If no more pair exists, single root may be returned
X and first order polynomial...
X """
X if len(self.conjs):
X r=self.conjs.pop(0)
X c=scipy.real(scipy.poly(r))
X return sos_n_or_d(r,c)
X elif len(self.real_pairs):
X r=self.real_pairs.pop(0)
X c=scipy.poly(r)
X c[1]=0 #Avoid non-zero c[1], possible due to numerical errors...
X return sos_n_or_d(r,c)
X elif len(self.reals)>1:
X r1=self.reals.pop(0)
X r2=self.reals.pop(0)
X c=scipy.poly((r1,r2))
X return sos_n_or_d(r,c)
X elif len(self.reals)==1:
X r=(self.reals.pop(0),)
X c=scipy.poly(r)
X return sos_n_or_d(r,c)
X else:
X r=None
X c=scipy.poly(1)
X return sos_n_or_d(r,c)
X
def rp_cmp(x,y):
X """
X This function is used for sorting of roots
X We sort the first - considering their real parts, then the imaginary
X """
X x=complex(x); y=complex(y) #Make sure, that the numbers are complex
X if x.real < y.real:
X return -1
X if x.real > y.real:
X return 1
X #Now compare the absolute value of the imaginary parts
X if abs(x.imag) < abs(y.imag):
X return -1
X if abs(x.imag) > abs(y.imag):
X return 1
X return 0
X
def pair_roots(r,tol=1000*scipy.MachAr().eps):
X """
X This function tries to pair roots
X """
X result = paired_roots()
X #First we go through the list, removing all real roots
X i=0
X rr=[] #list of real roots
X while i<len(r):
X if abs(scipy.imag(r[i]))<tol:
X rr.append(scipy.real(r.pop(i)))
X else:
X i+=1
X rr.sort(lambda x,y: cmp(abs(x),abs(y)))
X while len(rr)>0:
X print "real:"+str(rr[0])+" l:"+str(len(rr))
X #Check for the pair, considering the fact, that we may have multiple roots
X #Therefore we have to check all roots differing less than "tol"
X j=1
X while True:
X if ( j>=len(rr) ) or ( abs(rr[j])-abs(rr[0])>tol ) :
X #All roots with "almost the same" abs value checked!
X #No opposite root found. Add it as single root
X result.add_real([rr.pop(0),])
X break
X if abs(rr[j]+rr[0])<tol:
X #This is the opposite root, create the pair!
X b1 = rr.pop(j)
X b2 = rr.pop(0)
X result.add_real_pair(b1,b2)
X else:
X #Check next root
X j+=1
X #We sort the remainin roots in order of increasing real part, and then in order
X #of increasing imaginary part
X r.sort(rp_cmp)
X #After sorting we should be able to pair all roots as conjugate pairs
X #Unfortunately we can't assume that the neighbouring roots are conjugate.
X #Let's assume, that we have dual conjugate roots: a+bj, a+bj, a-bj, a-bj.
X #In this case due to numerical errors they may be ordered in different
X #ways...
X #Therefore the algorithm below is rather complicated...
X while len(r)>0:
X #We always work with the root r[0], as we remove paired roots from the list
X #print "conj:"+str(r[0])+" l:"+str(len(r))
X #Now we should find the conjugate
X j=1
X while True :
X if ( j>=len(r) ) or \
X ( abs(scipy.real(r[j]-r[0])) > tol ) or \
X (abs(scipy.imag(r[j]))-abs(scipy.imag(r[0])) > tol):
X #No conjugate found! It must be an error!!!
X #print "conj:"+str(r[0])+" l:"+str(len(r))
X raise "No matching conjugate!"
X else:
X if abs(r[j]-scipy.conj(r[0])) < tol:
X # This is really a conjugate root
X b1 = r.pop(j)
X b2 = r.pop(0)
X result.add_conjs(b1,b2)
X break
X else:
X #Check the next root
X j+=1
X return result
X
def iir2casc(tol,z_or_b, p_or_a, k=0):
X """
X This function converts the IIR filter described with the list of
X coefficients (b and a) or list with zeros, list of poles, and amplification k,
X into the cascaded form.
X The function returns the amplification coefficient k, and the list of
X second order sections...
X Currently no section zeros&poles grouping and section ordering
X is perfomred to optimize noise properties.
X """
X sections=[]
X # If k=0, it means, that k was not provided (k=0 makes no sense)
X # so we assume, that the input parameters are the filter's
X # b and a coefficients
X if k==0:
X #Convert the coefficients to the zpk form
X b=z_or_b
X a=p_or_a
X k=b[0]; b=b/b[0] #Is it right? What about the filters with b0=0?
X k/=a[0]; a=a/a[0]
X z=scipy.roots(b).tolist()
X p=scipy.roots(a).tolist()
X else:
X z=z_or_b
X p=p_or_a
X # OK. So here we have the zeroes, poles, and k
X # From this point the procedure is uniform
X # First we pair the complex zeroes and poles
X pz = pair_roots(z,tol)
X pp = pair_roots(p,tol)
X # Now we should group them and build the transfer functions...
X #
X # Here we group pairs, creating the second order sections.
X #We should replace it later with the more
X # inteligent grouping...
X while True:
X n=pz.get_sos_n_or_d()
X d=pz.get_sos_n_or_d()
X if not ( n.r or d.r ):
X # No more poles and zeroes!
X break
X sections.append(sos(n,d))
X return (k,sections)
X
# Below are the test routines
#
# The difficult case to test our roots pairing routine...
def test_roots_pairing():
X s=scipy.poly1d([1,1,-1,-2+1j,-2-1j,-2+1j,-2-1j,-2+2j,-2-2j],1)
X t=pair_roots(s.roots.tolist(),0.00001)
X return t
X
def test_iir2casc():
X (b,a)=scipy.signal.iirdesign([0.2, 0.5],[0.1, 0.6],0.1,40, 'ellip')
X print b
X print a
X casc=iir2casc(0.00001,b,a)
X print("k="+str(casc[0]))
X for s in casc[1]:
X print(s)
if __name__=='__main__':
X test_iir2casc()
SHAR_EOF
(set 20 09 05 09 15 30 25 'iir2casc.py'; eval "$shar_touch") &&
chmod 0600 'iir2casc.py'
if test $? -ne 0
then ${echo} 'restore of iir2casc.py failed'
fi
if ${md5check}
then (
${MD5SUM} -c >/dev/null 2>&1 || ${echo} 'iir2casc.py: MD5 check failed'
) << SHAR_EOF
802661fb2114afa736078422eca08c87 iir2casc.py
SHAR_EOF
else
test `LC_ALL=C wc -c < 'iir2casc.py'` -ne 8460 && \
${echo} 'restoration warning: size of iir2casc.py is not 8460'
fi
fi
if rm -fr ${lock_dir}
then ${echo} 'x - removed lock directory `'${lock_dir}\''.'
else ${echo} 'x - failed to remove lock directory `'${lock_dir}\''.'
exit 1
fi
exit 0