KG_is_back wrote: martinvicanek wrote:You can always factorize them (offline) into biquads, though.
I was always wondering how that is done.
Consider the Formant Filter from musicdsp.org (ported directly to FS):
Code: Select all
// Formant Filter
// http://www.musicdsp.org/showone.php?id=110
streamin in;
streamout out;
float out1, out2, out3, out4, out5, out6, out7, out8, out9, out10;
// vowel "A" coefficients:
float b0 = 8.11044e-06,
a1 = -8.943665402, a2 = 36.83889529, a3 = -92.01697887, a4 = 154.337906, a5 = -181.6233289,
a6 = 151.8651235, a7 = -89.09614114, a8 = 35.10298511, a9 = -8.388101016, a10 = 0.923313471;
out = b0*in - a1*out1 - a2*out2 - a3*out3 - a4*out4 - a5*out5 -
a6*out6 - a7*out7 - a8*out8 - a9*out9 - a10*out10;
out10 = out9;
out9 = out8;
out8 = out7;
out7 = out6;
out6 = out5;
out5 = out4;
out4 = out3;
out3 = out2;
out2 = out1;
out1 = out;
It is a 10th order allpole filter. Its direct implementation in FS' single precision is unstable. However, we can fix it by breaking it up in five 2nd order sections.
The transfer function of the formant filter is:
Code: Select all
b0
H(z) = --------------------------------------- (1)
1 + a1*z^-1 + a2*z^-2 + ... + a10*z^-10
We want to write it as:
Code: Select all
b0
H(z) = ---------------------------------------, (2)
(1 - z1/z)*(1 + z2/z)* ... *(1 - z10/z)
where the z1, z2, etc. are the poles of the filter. They are the roots of the denominator polynomial in equation (1). They can be real or complex conjugate pairs. There is a very useful root finder at
http://www.akiti.ca/PolyRootRe.html from which we get the following result:

- zeroes.png (61.83 KiB) Viewed 35450 times
Observe that the ten roots come in five complex conjugate pairs:
in equation (2) each such pair yields a section with real coefficients p and q:
Code: Select all
(1 - z1/z)*(1 - z2/z) = 1 + p*z^-1 + q*z^-2
with
The transfer function can be written as a product of 2nd order sections:
where
Code: Select all
b
H1(z) = ------------------- etc.
1 + p*z^-1 + q*z^-2
All that remains is to implement the 2nd order sections:
Code: Select all
// 2nd Order Section
streamin in;
streamin b;
streamin p;
streamin q;
streamout out;
float out1, out2;
out = b*in - p*out1 - q*out2;
out2 = out1;
out1 = out;
Below is a schematic for demonstration.