Frequency response of a complex IIR filter (Piotr Holonowicz )


Subject: Frequency response of a complex IIR filter
From:    Piotr Holonowicz  <piotr.holonowicz@xxxxxxxx>
Date:    Thu, 26 Jan 2012 15:02:41 +0100
List-Archive:<http://lists.mcgill.ca/scripts/wa.exe?LIST=AUDITORY>

This is a multi-part message in MIME format. --------------000809000204020909010803 Content-Type: text/html; charset=ISO-8859-1 Content-Transfer-Encoding: 7bit <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"> <html style="direction: ltr;"> <head> <meta content="text/html; charset=ISO-8859-1" http-equiv="Content-Type"> <title></title> <style>body p { margin-bottom: 0cm; margin-top: 0pt; } </style> <style>body p { margin-bottom: 0cm; margin-top: 0pt; } </style> </head> <body style="direction: ltr;" bidimailui-detected-decoding-type="latin-charset" bgcolor="#ffffff" text="#000000"> <p><br> </p> <meta content="text/html; charset=ISO-8859-1" http-equiv="Content-Type"> <title></title> <style>body p { margin-bottom: 0cm; margin-top: 0pt; }</style>Hi all!<br> I am trying to make a digital implementation of the One Zero Gammatone Filterbank (OZGF), however, I need the output to be complex instead of real. I have converted the filter equations by R. Lyon with the partial fraction expansion obtaining two conjugate poles, then treated it with the bilinear transform. The final stage is the coupled form of a complex IIR. In order to check the implementation, I would like to plot the frequency response of a single filter channel (= a cascade of three low-pass biquads + one bandpass biquad). But here is where I got a bit puzzled - for a real domain filter,&nbsp; the way to do it is to plot the magnitude and the phase of the FFT of the impulse response. How to do it however, with a complex output u+jv (&lt;=&gt; r exp (jw0)) ? My intuition says it might be done by computing y= r exp (w0) and then computing abs(fft(y)) and angle(fft(y)). Is this the correct way to do it or shall I treat the real and the imaginary output separately - in that case how to interpret the outcome? <br> I enclose the python prototype code, it requires only the numpy/scipy and the matplotlib(pylab) to run.<br> I'll be very thankful for help.<br> Best<br> Piotr Holonowicz<br> <br> PhD candidate at Music Technology Group<br> Universitat Pompeu Fabra<br> Barcelona,Spain<br> <br> </body> </html> --------------000809000204020909010803 Content-Type: text/x-python; name="bpbiquad.py" Content-Transfer-Encoding: 7bit Content-Disposition: attachment; filename="bpbiquad.py" #!/usr/bin/env python #****************************************************************************** # # Name: bpbigquad.py # Purpose: Bandpass biquad test. # Date: %date% at %time% # # # by Piotr Holonowicz <pholonow@xxxxxxxx> # Music Technology Group, University Pompeu Fabra, Barcelona, Spain # www.mtg.upf.es # ***************************************************************************** #imports import sys,os,glob from optparse import OptionParser import numpy import scipy from scipy import signal import pylab def bpbiquad(X, w0, Q, bp=False): """ Complex bandpass biquad IIR - plots the impulse response and the magnitude and phase Coupled form, obtained by the bilinear transform from its Laplace transfer function. w0 - omega0 - center frequency Q - quality factor """ d=1.0/Q #easier to compute Y = numpy.zeros(X.shape) mag = numpy.zeros(X.shape[0]) phs = numpy.zeros(X.shape[0]) om0 = numpy.tan(w0*0.5) s0 = numpy.sqrt(abs(d*d-4.0)) #bandpass or lowpass? if(bp): reZ0 = 0.5*om0 imZ0 = (om0 / d) * s0 #one zero else: reZ0 = 1.0 imZ0 = 1.0 #all pole reP0 = 0.5 * (om0 - d) imP0 = 0.5 * (om0 * s0) for i in range(X.shape[0]-1): Y[i, 0] = reZ0 * (reP0 * Y[i-1, 0] - imP0*Y[i-1, 1] + X[i-1, 0]) #last value is 0 , wrapping Y[i, 1] = imZ0 * (reP0 * Y[i-1, 1] - imP0*Y[i-1, 0] + X[i-1, 1]) mag[i] = numpy.sqrt(Y[i, 0] * Y[i, 0] + Y[i, 1] * Y[i, 1]) if(Y[i, 0]!=0.0): phs[i] = numpy.arctan(Y[i, 1] / Y[i, 0]) else: phs[i] = 0.0 return Y, mag, phs def magn(x): n = len(x)/2 return abs(scipy.fft(x))[:n] def phse(x): n = len(x)/2 return scipy.angle(scipy.fft(x))[:n] #main if __name__ == "__main__": x = numpy.zeros(128); x[0]=1.0; #impulse #x = scipy.signal.chirp(scipy.arange(100), 20.0, 1000.0, 20000.0) #make a complex version of the input signal, X[n] = [Re{x(n)}, Im{x(n)}] X = numpy.array([x, numpy.zeros(len(x))]).T w0 = scipy.pi/4.0; Q =2.5; Y, mag, phs = bpbiquad(X, w0, Q, bp=False); Y, mag, phs = bpbiquad(Y, w0, Q, bp=False); Y, mag, phs = bpbiquad(Y, w0, Q, bp=False); Y, mag, phs = bpbiquad(Y, w0, Q, bp=True); cplx = mag*scipy.exp(phs); pylab.close('all') xt = numpy.linspace(0.0, numpy.pi, len(x)) pylab.subplot(121) pylab.plot(mag) pylab.title("mag"); pylab.subplot(122) pylab.plot(phs) pylab.title("phase") #pylab.subplot(323) #pylab.plot(magn(cplx)) #pylab.title("Re_fftMag") #pylab.subplot(324) #pylab.plot(phse(cplx)) #pylab.title("Re_fftPhase") #pylab.subplot(325) #pylab.plot(magn(Y[:, 1])) #pylab.title("Im_fftMag") #pylab.subplot(326) #pylab.plot(phse(Y[:, 1])) #pylab.title("Im_fftPhase") pylab.show() --------------000809000204020909010803--


This message came from the mail archive
/var/www/postings/2012/
maintained by:
DAn Ellis <dpwe@ee.columbia.edu>
Electrical Engineering Dept., Columbia University