Hi all, I've been searching for the best way to do this programmically, but I can't seem to come up with a simple solution. I have a program that calculates filter coefficients correctly based on this z-transform: H(z) = (b0+b1*z^-1+b2*z^-2)/(1+a1*z^-1+a2*z^-2) I want to be able to display the magnitude response of this filter from 20Hz-20kHz at specified frequencies (48pts per octave) which I already have defined. Basically, I want to be able to pass in a frequency, lets say 250Hz, and be able to return a value based on what the filter is doing magnitude wise at that frequency assuming the sampling rate is 48kHz. Any help would be appreciated. Thanks, Dave

# Calculate frequency response of Z transform with known frequency bins

Started by ●March 6, 2007

Reply by ●March 6, 20072007-03-06

davidrhunt@gmail.com wrote:> Hi all, > > I've been searching for the best way to do this programmically, but I > can't seem to come up with a simple solution. > > I have a program that calculates filter coefficients correctly based > on this z-transform: > > H(z) = (b0+b1*z^-1+b2*z^-2)/(1+a1*z^-1+a2*z^-2) > > I want to be able to display the magnitude response of this filter > from 20Hz-20kHz at specified frequencies (48pts per octave) which I > already have defined. Basically, I want to be able to pass in a > frequency, lets say 250Hz, and be able to return a value based on what > the filter is doing magnitude wise at that frequency assuming the > sampling rate is 48kHz. > > Any help would be appreciated. > > Thanks, > Dave >Do you want to do this strictly inside your sampled-time world, or do you want to sample a real world signal, run it through the above filter, then reconstruct it back into a real world signal? For the strictly sampled time case you just set z = e^(1j*2*pi*250Hz/48kHz), then evaluate the amplitude of H(z). For the real world to real world case you need to take the anti-alias and reconstruction filters into account. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com Posting from Google? See http://cfaj.freeshell.org/google/ "Applied Control Theory for Embedded Systems" came out in April. See details at http://www.wescottdesign.com/actfes/actfes.html

Reply by ●March 6, 20072007-03-06

On Mar 6, 5:24 pm, davidrh...@gmail.com wrote:> Hi all, > > I've been searching for the best way to do this programmically, but I > can't seem to come up with a simple solution. > > I have a program that calculates filter coefficients correctly based > on this z-transform: > > H(z) = (b0+b1*z^-1+b2*z^-2)/(1+a1*z^-1+a2*z^-2) > > I want to be able to display the magnitude response of this filter > from 20Hz-20kHz at specified frequencies (48pts per octave) which I > already have defined. Basically, I want to be able to pass in a > frequency, lets say 250Hz, and be able to return a value based on what > the filter is doing magnitude wise at that frequency assuming the > sampling rate is 48kHz. > > Any help would be appreciated.2 years ago this came up. there is a numerical issue (having to do with cos(w) being so close to 1 when your frequencies are small) if you don't do this right. check out this thread, particularly posts 2 and 9: http://groups.google.com/group/comp.dsp/browse_frm/thread/8c0fa8d396aeb444/> Thanks,FWIW r b-j

Reply by ●March 7, 20072007-03-07

robert bristow-johnson wrote:> On Mar 6, 5:24 pm, davidrh...@gmail.com wrote: > 2 years ago this came up. there is a numerical issue (having to do > with cos(w) being so close to 1 when your frequencies are small) if > you don't do this right. check out this thread, particularly posts 2 > and 9: > > http://groups.google.com/group/comp.dsp/browse_frm/thread/8c0fa8d396aeb444/There's some nice info there in that thread. I must admit though that for filters where the critical frequency is very much lower than the sample rate, my preferred method of calculating frequency response is to take the FFT of the impulse response. This tends to capture rounding errors in calculation of the coefficients and/or in the topology of the filter much better than other methods and is less likely to give errors due to rounding during the freq. response calculation. Erik -- +-----------------------------------------------------------+ Erik de Castro Lopo +-----------------------------------------------------------+ "I'd rather not work with people who aren't careful. It's darwinism in software development." -- Linus Torvalds on the linux-kernel list

Reply by ●March 7, 20072007-03-07

On Mar 6, 11:22 pm, Erik de Castro Lopo <e...@mega-nerd.com> wrote:> > I must admit though that for filters where the critical frequency is > very much lower than the sample rate, my preferred method of calculating > frequency response is to take the FFT of the impulse response. This tends > to capture rounding errors in calculation of the coefficients and/or in the > topology of the filter much better than other methods and is less likely > to give errors due to rounding during the freq. response calculation.not that it makes any difference to me how you choose to do it, Erik, if it were me, i would prefer to be consistent in operation independent of where f0 goes (although we know there are numerical nasties when f0 << Fs). if you divide a0 into a1, a2, b0, b1, b2 and round to the precision used in your filter (i.e. a1, a2, b0, b1, b2 are exactly what your filter is using) the formula given should still be correct. but you *will* have to use a higher (double) precision in the frequency response code than is used in your filter so there are no overflows when f << Fs whether or not f0 << Fs. i think even this implementation of the frequency response curve is still far simpler than 1. evaluate impulse response in filter simulation 2. FFT impulse response, and 3. interpolate result of FFT to exponentially spaced frequencies for display. with the frequency response formula given, you can evaluate it at the discrete frequencies of interest (one for each pixel that you are plotting) whether they are linearly spaced in frequency or logarithmically spaced. r b-j

Reply by ●March 8, 20072007-03-08

robert bristow-johnson wrote:> not that it makes any difference to me how you choose to do it, Erik,It might not make any difference to me, but I might learn something :-).> if it were me, i would prefer to be consistent in operation > independent of where f0 goes (although we know there are numerical > nasties when f0 << Fs).Particularly for higher order filters. For really high order (36th order implemented in direct form in double floating point) f0 = 1/8 Fs already starts having problems while the impulse response method just works. Erik -- +-----------------------------------------------------------+ Erik de Castro Lopo +-----------------------------------------------------------+ "If dolphins are so smart, why do they live in igloos?" -Eric Cartman

Reply by ●March 8, 20072007-03-08

On Mar 8, 4:27 am, Erik de Castro Lopo <e...@mega-nerd.com> wrote:> robert bristow-johnson wrote: > > not that it makes any difference to me how you choose to do it, Erik, > > It might not make any difference to me, but I might learn something :-). > > > if it were me, i would prefer to be consistent in operation > > independent of where f0 goes (although we know there are numerical > > nasties when f0 << Fs). > > Particularly for higher order filters. For really high order (36th order > implemented in direct form in double floating point) f0 = 1/8 Fs already > starts having problems while the impulse response method just works.but how do you evaluate the impulse response? is that not just as nasty (even more so), numerically, as evaluating the terms of 20*log10[|H(e^jw)|] = 10*log10[ phi*( b0*b2*phi - (b1*(b0+b2) + 4*b0*b2) ) + (b0+b1+b2)^2 ] - 10*log10[ phi*( a0*a2*phi - (a1*(a0+a2) + 4*a0*a2) ) + (a0+a1+a2)^2 ] where phi = (2*sin(w/2))^2 and adding them up? the only numerical nasties happen when the terms with phi are so much smaller than the terms that they are added to (and that can be mitigated only with more bits in the words). when f0<<Fs nasty stuff happens in just calculating the coefficients (everybody suffers that problem). and more nasties happen in the implementation of the filter (very large quantities of opposite sign and nearly equal magnitude adding together to get a much smaller quantity and losing bits of precision as a result). this latter nasty is avoided in just evaluating the expression above for each biquad stage. and no FFT. that's my spin.> "If dolphins are so smart, why do they live in igloos?" -Eric Cartmanlast night i saw the latest S.P. with a stab at "Wheel of Fortune". the category was "People who annoy you" and the word was NAGGERS (but the "A" wasn't visible yet). wonderfully cynical politically incorrect stuff. South Park and The Daily Show are the only reasons i watch cable at all. r b-j

Reply by ●June 5, 20082008-06-05

Dear Sir I'm a master student in the academy of higher studies in Libya in the stage of writing the theses in the field of applied mathematics ,specificly (z-transform in Geophysics ), to be more specific I'm working on How to use and implement z-transform in oil industry ,But I have a problem with having previous studies, related web sites and more information to complete my these , the reason why I'm asking your help because we dont have enough books in this field so , I need your help to send me what could help me to complete my theses Looking forward to here from you Thanks alot omnia rabee