Calculate frequency response of Z transform with known frequency bins

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