DSPRelated.com
Forums

IIR Coefficient Calculation

Started by ger_lough April 17, 2014
On 4/20/14 7:21 PM, gyansorova@gmail.com wrote:
> > your example is fine since you sample 8 times the corner frequency. But you need to remind them that if you have a sampling frequency of say 3 times then you must use frequency pre-warping.
dunno exactly whom you're addressing. it it was me, the cookbook *does* account for frequency warping, both for location of the resonant (or "significant") frequency and for bandwidth. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
On Monday, April 21, 2014 12:14:34 PM UTC+12, robert bristow-johnson wrote:
> On 4/20/14 7:21 PM, gyansorova@gmail.com wrote: > > > > > > your example is fine since you sample 8 times the corner frequency. But you need to remind them that if you have a sampling frequency of say 3 times then you must use frequency pre-warping. > > > > dunno exactly whom you're addressing. it it was me, the cookbook *does* > > account for frequency warping, both for location of the resonant (or > > "significant") frequency and for bandwidth. > > > > -- > > > > r b-j rbj@audioimagination.com > > > > "Imagination is more important than knowledge."
I am looking at this link http://www.claysturner.com/dsp/Butterworth%20Filter%20Formulae.pdf which uses the Bilinear Transform. That is only accurate for larger sampling rates wrt the cut-off. You must be refering to another link. Also I have not seen the cot(pi*fc/fs) bit before - is that what I am missing? I use 2/T or T/2 depending on which way I am converting. Do you have a proof for the cot() bit?
On 4/20/14 8:54 PM, gyansorova@gmail.com wrote:
> > > I am looking at this link > http://www.claysturner.com/dsp/Butterworth%20Filter%20Formulae.pdf > > which uses the Bilinear Transform. That is only accurate for larger sampling rates wrt the cut-off. You must be referring to another link.
i was looking for an "example", and i didn't see that in Clay's paper.
> Also I have not seen the cot(pi*fc/fs) bit before - is that what I am missing? I use 2/T or T/2 depending on which way I am converting.
i think it's normalized out. when you do the bilinear transform *with* compensation for frequency warping, and apply it to analog prototypes that have *normalized* resonant frequency, i think the substitution s <-- cot(pi*fc/fs) * (1 - z^-1)/(1 + z^-1) is sufficient.
> Do you have a proof for the cot() bit?
i am not Clay. but in my experience, this cotangent comes directly from applying frequency warping to the "fc" frequency before doing the bilinear transform substitution. in *my* treatment of this, i immediately get rid of the cotangent (and cotangent squared) with some trig substitutions that are spelled out in the cookbook. but, as best as i can tell, Clay *has* dealt with the frequency warping issue in relocating the significant frequency. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
On Monday, April 21, 2014 1:26:41 PM UTC+12, robert bristow-johnson wrote:
> On 4/20/14 8:54 PM, gyansorova@gmail.com wrote: > > > > > > > I am looking at this link > > > http://www.claysturner.com/dsp/Butterworth%20Filter%20Formulae.pdf > > > > > > which uses the Bilinear Transform. That is only accurate for larger sampling rates wrt the cut-off. You must be referring to another link. > > > > i was looking for an "example", and i didn't see that in Clay's paper. > > > > > Also I have not seen the cot(pi*fc/fs) bit before - is that what I am missing? I use 2/T or T/2 depending on which way I am converting. > > > > i think it's normalized out. when you do the bilinear transform *with* > > compensation for frequency warping, and apply it to analog prototypes > > that have *normalized* resonant frequency, i think the substitution > > > > s <-- cot(pi*fc/fs) * (1 - z^-1)/(1 + z^-1) > > > > is sufficient. > > > > > Do you have a proof for the cot() bit? > > > > i am not Clay. but in my experience, this cotangent comes directly from > > applying frequency warping to the "fc" frequency before doing the > > bilinear transform substitution. in *my* treatment of this, i > > immediately get rid of the cotangent (and cotangent squared) with some > > trig substitutions that are spelled out in the cookbook. > > > > but, as best as i can tell, Clay *has* dealt with the frequency warping > > issue in relocating the significant frequency. > > > > -- > > > > r b-j rbj@audioimagination.com > > > > "Imagination is more important than knowledge."
If you use normalised polynomial you must first do a low pass to (say) low pass transform before you convert to digital. at least that's the method I know. eg s->s/wc or s->wc/s for low pass to high pass. the normalised prototype is always low pass.
On Friday, April 18, 2014 8:32:13 AM UTC+12, ger_lough wrote:
> I am having difficulty finding any info on how to calculate digital IIR > > filter coefficients. I know this can be done in MATLAB, for example, but I > > am interested in changing IIR specs on the fly. I have done this before > > using the formulae provided by Robert Bristow-Johnson on musicdsp.org, but > > these equations only cover a 2nd order biquad. So, for example, how can I > > calculate the coefficients for a Butterworth (low-pass, band-pass, > > high-pass, etc) of order 32 using the same input as provided by RBJ on > > musicdsp.org (Fs, f0 and Q). I need formulae so that the coefficients can > > be updated on the fly. Also, it would be much appreciated if somebody can > > point me to some literature specifically on this subject (i.e. how to > > calculate iir coefficients for various IIR types). > > > > > > _____________________________ > > Posted through www.DSPRelated.com
Just to mention the way I know from the textbook. I assume you know the order of the filter and for simplicity make it a low pass. You need the ripple in the passband - to find the ripple factor eps which becomes unity in the 3dB case which is usually quoted. Otherwise the poles of the analogue filter have radius (1/eps)^1/n where eps is ripple factor and n is order. Knowing the radius of the poles and the phase angles of the poles you can then calculate the normalised low-pass prototype. You get the classic Butterworth polynomials in the 3dB case (as in the case of the notes that were posted) but not otherwise. Then you transform low pass to low pass or low pass to high pass, low pass to bandpass or low pass to bandstop using 1 of 4 transforms in the textbooks. Then you have your analogue Butterworth filter. For digital we use the Bilinear TF making sure to sample 8-10 times higher than the highest frequency otherwise you have to piss about with pre-warping.
On 4/21/14 2:00 AM, gyansorova@gmail.com wrote:
>> > If you use normalised polynomial you must first do a low pass to (say) low pass transform before you convert to digital. at least that's the method I know. > > eg s->s/wc or s->wc/s for low pass to high pass. the normalised prototype is always low pass.
and in all cases, you can integrate the T into that and also integrate the fudge factor we do to compensate for frequency warping. do all three in one operation since it's just scaling. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
On Tuesday, April 22, 2014 1:59:08 AM UTC+12, robert bristow-johnson wrote:
> On 4/21/14 2:00 AM, gyansorova@gmail.com wrote: > > >> > > > If you use normalised polynomial you must first do a low pass to (say) low pass transform before you convert to digital. at least that's the method I know. > > > > > > eg s->s/wc or s->wc/s for low pass to high pass. the normalised prototype is always low pass. > > > > > > and in all cases, you can integrate the T into that and also integrate > > the fudge factor we do to compensate for frequency warping. do all > > three in one operation since it's just scaling. > > > > -- > > > > r b-j rbj@audioimagination.com > > > > "Imagination is more important than knowledge."
oh very good
On Sunday, April 20, 2014 8:54:54 PM UTC-4, gyans...@gmail.com wrote:
> On Monday, April 21, 2014 12:14:34 PM UTC+12, robert bristow-johnson wrote: > > > On 4/20/14 7:21 PM, gyansorova@gmail.com wrote: > > > > > > > > > > > > > > your example is fine since you sample 8 times the corner frequency. But you need to remind them that if you have a sampling frequency of say 3 times then you must use frequency pre-warping. > > > > > > > > > > > > dunno exactly whom you're addressing. it it was me, the cookbook *does* > > > > > > account for frequency warping, both for location of the resonant (or > > > > > > "significant") frequency and for bandwidth. > > > > > > > > > > > > -- > > > > > > > > > > > > r b-j rbj@audioimagination.com > > > > > > > > > > > > "Imagination is more important than knowledge." > > > > I am looking at this link > > http://www.claysturner.com/dsp/Butterworth%20Filter%20Formulae.pdf > > > > which uses the Bilinear Transform. That is only accurate for larger sampling rates wrt the cut-off. You must be refering to another link. > > > > Also I have not seen the cot(pi*fc/fs) bit before - is that what I am missing? I use 2/T or T/2 depending on which way I am converting. Do you have a proof for the cot() bit?
Actually the cutoff freq is correct. That is in the cot(pi*f) part of the transform. The bilinear transform maps three points between the Laplace and Z domains. Namely 0 -> 0, infty -> half of the sampling rate, and finally 1 rad/sec -> to your digital cutoff freq. That is why I gave a formula for the Butterworths normalized to 1 rad/sec. At the cutoff, you will have a magnitude response of sqrt(1/2) e.g., approx. -3dB. Clay
On Monday, April 21, 2014 2:00:41 AM UTC-4, gyans...@gmail.com wrote:
> On Monday, April 21, 2014 1:26:41 PM UTC+12, robert bristow-johnson wrote: > > > On 4/20/14 8:54 PM, gyansorova@gmail.com wrote: > > > > > > > > > > > > > > > I am looking at this link > > > > > > > http://www.claysturner.com/dsp/Butterworth%20Filter%20Formulae.pdf > > > > > > > > > > > > > > which uses the Bilinear Transform. That is only accurate for larger sampling rates wrt the cut-off. You must be referring to another link. > > > > > > > > > > > > i was looking for an "example", and i didn't see that in Clay's paper. > > > > > > > > > > > > > Also I have not seen the cot(pi*fc/fs) bit before - is that what I am missing? I use 2/T or T/2 depending on which way I am converting. > > > > > > > > > > > > i think it's normalized out. when you do the bilinear transform *with* > > > > > > compensation for frequency warping, and apply it to analog prototypes > > > > > > that have *normalized* resonant frequency, i think the substitution > > > > > > > > > > > > s <-- cot(pi*fc/fs) * (1 - z^-1)/(1 + z^-1) > > > > > > > > > > > > is sufficient. > > > > > > > > > > > > > Do you have a proof for the cot() bit? > > > > > > > > > > > > i am not Clay. but in my experience, this cotangent comes directly from > > > > > > applying frequency warping to the "fc" frequency before doing the > > > > > > bilinear transform substitution. in *my* treatment of this, i > > > > > > immediately get rid of the cotangent (and cotangent squared) with some > > > > > > trig substitutions that are spelled out in the cookbook. > > > > > > > > > > > > but, as best as i can tell, Clay *has* dealt with the frequency warping > > > > > > issue in relocating the significant frequency. > > > > > > > > > > > > -- > > > > > > > > > > > > r b-j rbj@audioimagination.com > > > > > > > > > > > > "Imagination is more important than knowledge." > > > > If you use normalised polynomial you must first do a low pass to (say) low pass transform before you convert to digital. at least that's the method I know. > > > > eg s->s/wc or s->wc/s for low pass to high pass. the normalised prototype is always low pass.
You may transform the filter type (e.g. low pass to high pass by letting s = 1/s) in the Laplace domain and then apply the bilinear transform. Example low pass -> H(s) = 1/(1+s) transform by replacing s with 1/s H(1/s) = 1/(1+1/s) = s/(1+s) high pass with cutoff at 1 rad/sec. Chebyshev filters work exactly the same way. Instead of the poles being on the unit circle, they are instead on an ellipse whose eccentricity is related to the oscillation factor for the Chebyshev filter. I'll have to look up the formula if you need it. Any good text on analog filter theory should contain these formulae. Clay
On Thursday, April 17, 2014 4:32:13 PM UTC-4, ger_lough wrote:
> I am having difficulty finding any info on how to calculate digital IIR > > filter coefficients. I know this can be done in MATLAB, for example, but I > > am interested in changing IIR specs on the fly. I have done this before > > using the formulae provided by Robert Bristow-Johnson on musicdsp.org, but > > these equations only cover a 2nd order biquad. So, for example, how can I > > calculate the coefficients for a Butterworth (low-pass, band-pass, > > high-pass, etc) of order 32 using the same input as provided by RBJ on > > musicdsp.org (Fs, f0 and Q). I need formulae so that the coefficients can > > be updated on the fly. Also, it would be much appreciated if somebody can > > point me to some literature specifically on this subject (i.e. how to > > calculate iir coefficients for various IIR types). > > > > > > _____________________________ > > Posted through www.DSPRelated.com
Not sure if this is what you are looking for: A. G. Constantinides, "Spectral transformations for digital filters", Proc. Inst. Elect. Eng., vol. 117, no. 8, pp.1585 -- 1590, 1970. If you have DSP System Toolbox, this algorithm is implemented in dsp.VariableBandwidthIIRFilter System object: http://www.mathworks.com/help/dsp/ref/dsp.variablebandwidthiirfilter-class.html The parameters like PassbandFrequency/CenterFrequency/Bandwidth are tunable on the fly.