Hi everyone! Right now I am programing an application with Visual C++ for IIR Filters and everything is working fine, I mean the coefficients of the filters are correct, but I don't really know what I should do for plotting the magnitude VS frequency response. If I pretend to use my recursive formula (from bilinear transformation) I need a hell of calculations as the filter needs some time until it starts providing the right values, and I should store all this values for each frequency so I can determine the max point. And if I choose a high sample rate... I can always plot the analog transfer function for the filter instead of the digital as it is REALLY close, but it is not EXACT. It would be great if you have some ideas or suggestions!! Thanks a lot to everyone!!!
IIR Filters: How can I plot an EXACT Bode Diagram?
Started by ●September 20, 2008
Reply by ●September 20, 20082008-09-20
On 20 Sep, 14:11, "Calculadora" <pablo.martinez.lina...@gmail.com> wrote:> Hi everyone! > Right now I am programing an application with Visual C++ for IIR Filters > and everything is working fine, I mean the coefficients of the filters are > correct, but I don't really know what I should do for plotting the > magnitude VS frequency response. > If I pretend to use my recursive formula (from bilinear transformation) I > need a hell of calculations as the filter needs some time until it starts > providing the right values, and I should store all this values for each > frequency so I can determine the max point. And if I choose a high sample > rate... > I can always plot the analog transfer function for the filter instead of > the digital as it is REALLY close, but it is not EXACT. > It would be great if you have some ideas or suggestions!! > Thanks a lot to everyone!!!I assume you want to plot the frequency resonse of th efilter? If so, why not compute the frequency response directly from the coefficients? How I would do this in matlab: N=512; Hw = 20*log10(abs(fft(b,N)./fft(a,N))); Rune
Reply by ●September 20, 20082008-09-20
On Sep 20, 8:48 am, Rune Allnor <all...@tele.ntnu.no> wrote:> > I assume you want to plot the frequency resonse of th efilter? > If so, why not compute the frequency response directly from > the coefficients? > > How I would do this in matlab: > > N=512; > Hw = 20*log10(abs(fft(b,N)./fft(a,N))); >wasn't freqz() s'pose to do this also? i know what fft(a) does, but what does fft(a,N) do? does it zero-pad a to be length N and then FFT it? i s'pose i could just look it up. r b-j
Reply by ●September 20, 20082008-09-20
On 20 Sep, 16:45, robert bristow-johnson <r...@audioimagination.com> wrote:> On Sep 20, 8:48 am, Rune Allnor <all...@tele.ntnu.no> wrote: > > > > > I assume you want to plot the frequency resonse of th efilter? > > If so, why not compute the frequency response directly from > > the coefficients? > > > How I would do this in matlab: > > > N=512; > > Hw = 20*log10(abs(fft(b,N)./fft(a,N))); > > wasn't freqz() s'pose to do this also?Maybe. I prefer to roll my own instead of using canned matlab routines, partially because I would like to know exactly what is going on, but mostly because I don't want to pay $1000 for the SP toolbox which I won't use anyway.> i know what fft(a) does, but what does fft(a,N) do? �does it zero-pad > a to be length N and then FFT it?Yep. Rune
Reply by ●September 20, 20082008-09-20
On Sat, 20 Sep 2008 07:11:06 -0500, "Calculadora" <pablo.martinez.linares@gmail.com> wrote:>I mean the coefficients of the filters are >correct, but I don't really know what I should do for plotting the >magnitude VS frequency response.This was discussed right here about a month ago. You can compute the frequency response directly from the transfer function coefficients: z = exp(sT) T = 1/(sampling frequency) s = sigma + j�omega Sigma equals zero and omega equals (2PI�frequency) when computing frequency response. exp(j�omega�T) = cos(omega�T) + j�sin(omega�T) Plug z = cos(omega�T) + j�sin(omega�T) into your transfer function and solve for magnitude and phase at whatever frequencies you desire. -- Greg
Reply by ●September 20, 20082008-09-20
On Sat, 20 Sep 2008 07:45:07 -0700 (PDT), robert bristow-johnson <rbj@audioimagination.com> wrote:>i know what fft(a) does, but what does fft(a,N) do? does it zero-pad >a to be length N and then FFT it?IIRC, it truncates length(a) to N if length(a) > N, and zeropads if length(a) < N. Greg
Reply by ●September 20, 20082008-09-20
Greg Berchin wrote:> On Sat, 20 Sep 2008 07:11:06 -0500, "Calculadora" > <pablo.martinez.linares@gmail.com> wrote: > >> I mean the coefficients of the filters are >> correct, but I don't really know what I should do for plotting the >> magnitude VS frequency response. > > This was discussed right here about a month ago. You can compute the > frequency response directly from the transfer function coefficients: > > z = exp(sT) > > T = 1/(sampling frequency) > > s = sigma + j�omega > > Sigma equals zero and omega equals (2PI�frequency) when computing > frequency response. > > exp(j�omega�T) = cos(omega�T) + j�sin(omega�T) > > Plug z = cos(omega�T) + j�sin(omega�T) into your transfer function and > solve for magnitude and phase at whatever frequencies you desire. >If you're doing real-world signal processing and if your bandwidth is close to your sampling rates you'll probably also want to factor in the effect of the DAC and possibly the reconstruction filters. Barring aliasing effects this is easy -- just multiply the above-derived frequency response with the frequency response of the zero-order-hold model of the DAC, then with the reconstruction filter's response, if applicable. -- Tim Wescott Wescott Design Services http://www.wescottdesign.com Do you need to implement control loops in software? "Applied Control Theory for Embedded Systems" gives you just what it says. See details at http://www.wescottdesign.com/actfes/actfes.html
Reply by ●September 24, 20082008-09-24
>On Sat, 20 Sep 2008 07:11:06 -0500, "Calculadora" ><pablo.martinez.linares@gmail.com> wrote: > >>I mean the coefficients of the filters are >>correct, but I don't really know what I should do for plotting the >>magnitude VS frequency response. > >This was discussed right here about a month ago. You can compute the >frequency response directly from the transfer function coefficients: > >z = exp(sT) > >T = 1/(sampling frequency) > >s = sigma + j�omega > >Sigma equals zero and omega equals (2PI�frequency) when computing >frequency response. > >exp(j�omega�T) = cos(omega�T) + j�sin(omega�T) > >Plug z = cos(omega�T) + j�sin(omega�T) into your transfer functionand>solve for magnitude and phase at whatever frequencies you desire. > >-- Greg >Hi everyone and thanks for the replies!! But Greg, when you suggest to plug z in my transfer function, is it in the DC Gain?? I think I didn't made myself very clear... The problem I found is that you need some time untill the recursive Transfer function starts giving the proper values for the frequency desired, and then I would have to simulate a sine with the desired frequency and amplitude before I get the right value. Maybe my aproach is incorrect, so I will try to do this you suggest (despite I don't see it very clear :S). Thanks a lot to everyone!
Reply by ●September 24, 20082008-09-24
Calculadora wrote:> But Greg, when you suggest to plug z in my transfer function, is it in the > DC Gain??No, z is the z-transform operator. You indicated that you have an IIR system. In the discrete time domain, such a system is represented by a transfer function of the form: H(z) = [b0 + b1z^(-1) + ... + bnz^(-n)] / [1 + a1z^(-1) + ... + adz^(- d)] Plug in the values for z from my previous message and solve for H(z) at each frequency of interest. H(z) will be a complex number, with a real part and an imaginary part. Perform a rectangular to polar conversion, and you will have magnitude and phase.> I think I didn't made myself very clear... The problem I found is that you > need some time untill the recursive Transfer function starts giving the > proper values for the frequency desired, and then I would have to simulate > a sine with the desired frequency and amplitude before I get the right > value.It looks like you're implementing the time domain difference equation and applying a sinewave of the appropriate frequency, waiting for the startup transient to decay to negligible level, and observing the magnitude and phase of the output sinewave. That's a truly brute- force method that doesn't work very well. Greg
Reply by ●September 25, 20082008-09-25
>Calculadora wrote: > >> But Greg, when you suggest to plug z in my transfer function, is it inthe>> DC Gain?? > >No, z is the z-transform operator. > >You indicated that you have an IIR system. In the discrete time >domain, such a system is represented by a transfer function of the >form: > >H(z) = [b0 + b1z^(-1) + ... + bnz^(-n)] / [1 + a1z^(-1) + ... + adz^(- >d)] > >Plug in the values for z from my previous message and solve for H(z) >at each frequency of interest. H(z) will be a complex number, with a >real part and an imaginary part. Perform a rectangular to polar >conversion, and you will have magnitude and phase. > >> I think I didn't made myself very clear... The problem I found is thatyou>> need some time untill the recursive Transfer function starts givingthe>> proper values for the frequency desired, and then I would have tosimulate>> a sine with the desired frequency and amplitude before I get the right >> value. > >It looks like you're implementing the time domain difference equation >and applying a sinewave of the appropriate frequency, waiting for the >startup transient to decay to negligible level, and observing the >magnitude and phase of the output sinewave. That's a truly brute- >force method that doesn't work very well. > >Greg >All Right! As you said what I was trying was totally crazy. But it is clear now how I should do it. Thanks a lot!
