Confused by 1pole IIR code

[gordonnn]

I came across the code posted here:
http://www.musicdsp.org/en/latest/Filte ... lp-hp.html

here is the code in question, there are some revisions posted in the link as well:

Code: Select all

void SetLPF(float fCut, float fSampling)
{
    float w = 2.0 * fSampling;
    float Norm;

    fCut *= 2.0F * PI;
    Norm = 1.0 / (fCut + w);
    b1 = (w - fCut) * Norm;
    a0 = a1 = fCut * Norm;
}

void SetHPF(float fCut, float fSampling)
{
    float w = 2.0 * fSampling;
    float Norm;

    fCut *= 2.0F * PI;
    Norm = 1.0 / (fCut + w);
    a0 = w * Norm;
    a1 = -a0;
    b1 = (w - fCut) * Norm;
}

Where
out[n] = in[n]*a0 + in[n-1]*a1 + out[n-1]*b1;

I was interested because seeing the use of exp() omitted when setting the cutoff seemed like it might be faster, as I want to be able to modulate the cutoff despite it being a simple filter.

How are the coefficients being set here? Is it something like an approximation to using exp by using two terms of a Taylor series?

I was wondering, because it looks completely unfamiliar to me.

Thanks for any help

[mystran]

Since we put a zero at either z=1 (HP) or z=-1 (LP), these are BLT (bilinear transform) and hence the correct tuning formula would be trigonometric, with the most obvious formula using tan(pi*fCut/fSampling).

As it turns out, tan(x) is approximately equal to x for small values of x, such that for low-frequencies, you can skip the tan() without too much error. After you make this approximation, the above coefficient formulas can probably (sorry, I'm not in a mood to really check properly) be obtained from the BLT by algebraic manipulation, although I honestly can't remember if I've ever seen it written quite like that.

[gordonnn]

It's true that omitting tan seems to work fine.
I've seen a common formula for a one pole filter being: b = tan(pi*fcut/sampling) and a = 1.0 - b, for out[n] = b*in[n] + a*out[n-1]
It looks like the pasted code approximates the same formula for, but with 2pi instead of pi, and approximates setting a to 1-b.
It also sounds about the same as lop~ in puredata which is a one pole filter with the simple formula, but also using 2pi instead of pi, and also since in the pasted code a1 == a0, it's practically the same cutoff and same sound as lop~ just louder.
Experimenting in PureData with half the cutoff set in lop~ sounds exactly like using pi instead of 2pi for the common formula..

I'm not sure what to make of all that in terms of what is a correct way to set coefficients. Does any particular approach strike you as being more accurate, and does the zero term do anything noticeable besides boost the signal? Sorry if these are vague questions, I probably won't be able to get a sense of "correct" without having taken differential equations anyway so I can really understand the math fully.

[mystran]

It's true that omitting tan seems to work fine.

It works for low-frequencies, but for higher frequencies it under-estimates the cutoff. If you need accurate tuning higher up, you really need at least an approximation of tan().

I'm not sure what to make of all that in terms of what is a correct way to set coefficients. Does any particular approach strike you as being more accurate, and does the zero term do anything noticeable besides boost the signal? Sorry if these are vague questions, I probably won't be able to get a sense of "correct" without having taken differential equations anyway so I can really understand the math fully.

Your question is vague enough that in order to really answer it properly would basically take a full tutorial into filter design, but basically the transfer function of an analog filter is a rational function in (complex) s (which is the derivative operator in time-domain) in the Laplace-transform domain where the Fourier transform (of the filter impulse response) is found on the imaginary axis of the s-plane. Similarly the transfer function of a digital filter is a rational function in (complex) z^-1 (which is the unit-delay in discrete time-domain) in the z-transform domain, where the DFT is found on the unit circle.

In both cases, placing zeroes or poles at various points on the complex plane will attenuate or boost (respectively) those frequencies (ie. the points on the imaginary axis or unit circle) that are "close" to the location of the zeroes or poles.

Generally speaking, if you want to learn this stuff, I suggest learning about analog filters first. You really don't need to worry about electronic components, but working with s-plane transfer functions is a lot easier. The bilinear transform then is basically what you get when you numerically integrate the continuous-time differential equations using trapezoidal integration.

One can actually do such integration literally (which for historical reasons is usually called ZDF filter design around here, see "The Art of VA Filter Design"), or you can substitute the trapezoidal differentiator in the analog transfer function as s<-2/T*(z - 1)/(z+1) and then simplify the result into a rational in z^-1 (or realistically, let a computer algebra system simplify it for you; doing this stuff manually will drive you nuts pretty quickly).

The tan() tuning basically amounts to "pre-warping" the analog frequencies in such a way that the unit-frequency (which is what is usually assumed whenever there isn't any explicit tuning parameter) ends up (after the transform) at the desired digital frequency.

I realise this probably won't really answer your questions very well, but hopefully it will give you an idea of what is involved and what to search for when you want to learn more.

[matt42]

Though it gets rightly criticised for poor coding examples I feel that Designing Audio Effect Plugins in C++ is probably one of the best resource for a beginner to get a basic grasp of digital filter theory. Not perfect with a lot of unnecessary repetition, but the theory is reasonably well explained for those without extensive maths backgrounds.

[juha_p]

Code: Select all

float lpf_MZTapprox(float fs, float fc){

// fs - in Hz (tested only @  44100)
// fc - in µs, if Hz used then convert to µs with formula 1/(2*pi*fc) 
	float x = 1.0/(fs*fc);
	// a0 = 1.0;
	// float a1 = -(0.99999999999999999+0.99999999999062503*x+0.49999750000703124*x*x)
	float a1 = -(1 + x + 0.5*x*x);
	float b0 = 1.0 + a1;
	// b1 = 0.0; 
}

float lpf_MZT(float fs, float fc){
// fs - in Hz
// fc - in Hz
	//float a0 = 1.0;
	float a1 = -exp(-1.0/(fs*(1/(2 * pi * fc)))); 
	float b0 = 1.0-a1;
	//float b1 = 0.0; 
}

float lpf_BLT(float fs, float fc){
    float w0 = 2*M_PI*(fc/fs);
    float b0 =   sin(w0);
    float b1 =   sin(w0);
    float a0 =   cos(w0) + sin(w0) + 1.0;
    float a1 =   sin(w0) - cos(w0) - 1.0;
}

By my testing it's about 20 times faster than original exp() based version and the accuracy is quite good ... knowing what MZT based LPF can be.

Idea used in "lpf_MZTapprox" is to use approximation error to cancel the error in MZT method ... when higher degree polynomial is used the resulting filter closes the original MZT/IIM responses.
See the values for s, c and v, v2: https://www.desmos.com/calculator/ygwscpxu1r

Someone with better math skills can easily improve this idea by finding better polynomial with suitable error to cancel the error in MZT method.

EDIT: Added code for "BLT" and std MZT version (+ updated the plot).
EDIT2: Added some info of the idea behind lpf_MZTapprox -function.

[gordonnn]

Really helpful, thanks! I will put // credit: juha_p if I use this one
Hopefully I will be able to work through the steps going from a transfer function to filter coefficients eventually but this is very useful.