Forums

Calculate frequency response of Z transform with known frequency bins

Started by Unknown March 6, 2007
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

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
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
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
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
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
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 Cartman
last 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
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