Charles Henry wrote:
wow - this is a bit above my head at the moment. i can go from equations in C to pd, but only if i have the equations themselves to look at...
It's easy. I've done it half a dozen times. For example, http://lists.puredata.info/pipermail/pd-list/2007-01/046315.html
Look at the help file for biquad~ in Pd. The equations are there.
last time i looked at the equations in biquad~'s help patch, my brain melted.
Here's the transfer function: Y/X=(ff1 + ff2*z^(-1) + ff3*z^(-2))/(1 - fb1*z^(-1) - fb2*z^(-2))
or
Y/X=(ff1*z^2 + ff2*z^1 + ff3)/(1*z^2 - fb1*z^1 - fb2)
All that you need to do to find the complex zeros is to use quadratic formula on the numerator (ff1*z^2 + ff2*z^1 + ff3)
and do the same to find the complex poles on (1*z^2 - fb1*z^1 - fb2)
dude, that's not 'easy'. i can barely remember how the quadratic formula works with pen and paper, let alone in C, and let alone to the point where i could confidently transfer from C to Pd. i'd have no idea where the bugs were... and indeed, no idea
and once i've got these 'complex zeros' and 'complex poles' i don't have any idea how to plug them into czero~ and cpole~ (or is it rzero~ and rpole~?) objects.
On Mon, Sep 15, 2008 at 7:49 AM, Damian Stewart damian@frey.co.nz wrote:
dude, that's not 'easy'. i can barely remember how the quadratic formula works with pen and paper, let alone in C, and let alone to the point where i could confidently transfer from C to Pd. i'd have no idea where the bugs were... and indeed, no idea
Use pen and paper. Then, take your results and make implementations in Pd (or C, your choice). This is a good exercise, so I'll leave it to you to do the derivation. I'm more than happy to help. I do this stuff all the time.
- Solve for x in
a*x^2 + b*x + c =0 (use "completing the square") This is how you derive the quadratic formula.
- Then, solve
ff1*z^2 + ff2*z^1 + ff3 = 0 and 1*z^2 - fb1*z^1 - fb2 = 0 using the result from (1)
- Then, solve for fb1, fb2, ff1, ff2, ff3 in terms of the poles/zero
locations from (2)
Then, you've done all the work you need to freely exchange between both methods. Now, as for filter design, there's a bit more work to do (save that for later).
Chuck
Hallo, Charles Henry hat gesagt: // Charles Henry wrote:
On Mon, Sep 15, 2008 at 7:49 AM, Damian Stewart damian@frey.co.nz wrote:
dude, that's not 'easy'. i can barely remember how the quadratic formula works with pen and paper, let alone in C, and let alone to the point where i could confidently transfer from C to Pd. i'd have no idea where the bugs were... and indeed, no idea
Use pen and paper. Then, take your results and make implementations in Pd (or C, your choice). This is a good exercise, so I'll leave it to you to do the derivation. I'm more than happy to help. I do this stuff all the time.
- Solve for x in
a*x^2 + b*x + c =0 (use "completing the square") This is how you derive the quadratic formula.
Okay, that's easy (or can be looked up e.g. here: http://www.purplemath.com/modules/sqrquad2.htm )
- Then, solve
ff1*z^2 + ff2*z^1 + ff3 = 0 and 1*z^2 - fb1*z^1 - fb2 = 0 using the result from (1)
Easy now as well.
- Then, solve for fb1, fb2, ff1, ff2, ff3 in terms of the poles/zero
locations from (2)
After (2) I have something like this:
x = (-b +- sqrt(b*b - 4*a*c)) / 2*a
zf = (-ff2 +- sqrt(ff2*ff2 - 4*ff1*ff3)) / 2*ff1 zb = (-fb1 +- sqrt(fb1*fb1 - 4*fb2)) / 2
This gives me four complex numbers zf+, zf-, zb+ and zb-, right?
Now in a biquad~-clone with czero~/cpole~, the values for fb1, fb2, ff1, ff2, ff3 would be known so I also could calculate zf and zb now.
But what I still don't get is how to feed this into a patch as below as coefficients to make a biquad~ clone:
[*~] ? | | ? ? [czero~] | | | | ? ? [czero~] | | | | ? ? [cpole~] | | | | ? ? [cpole~]
This still has several question marks. :(
Ciao
Frank
Hallo, Frank Barknecht hat gesagt: // Frank Barknecht wrote:
1*z^2 - fb1*z^1 - fb2 = 0
zb = (-fb1 +- sqrt(fb1*fb1 - 4*fb2)) / 2
Sorry, the one above is wrong, of course. Given the minus-signs in the original transfer formula it should be:
zb = (fb1 +- sqrt(fb1*fb1 + 4*fb2)) / 2
if I'm right.
Ciao
Frank
Frank Barknecht wrote:
Hallo, Frank Barknecht hat gesagt: // Frank Barknecht wrote:
1*z^2 - fb1*z^1 - fb2 = 0
zb = (-fb1 +- sqrt(fb1*fb1 - 4*fb2)) / 2
Sorry, the one above is wrong, of course. Given the minus-signs in the original transfer formula it should be:
zb = (fb1 +- sqrt(fb1*fb1 + 4*fb2)) / 2
if I'm right.
see, Chuck, if Frank gets it wrong what chance do i have? ;-)
d
damian stewart | skype: damiansnz | damian@frey.co.nz frey | live art with machines | http://www.frey.co.nz
Hallo, Damian Stewart hat gesagt: // Damian Stewart wrote:
Frank Barknecht wrote:
Hallo, Frank Barknecht hat gesagt: // Frank Barknecht wrote:
1*z^2 - fb1*z^1 - fb2 = 0
zb = (-fb1 +- sqrt(fb1*fb1 - 4*fb2)) / 2
Sorry, the one above is wrong, of course. Given the minus-signs in the original transfer formula it should be:
zb = (fb1 +- sqrt(fb1*fb1 + 4*fb2)) / 2
if I'm right.
see, Chuck, if Frank gets it wrong what chance do i have? ;-)
;)
Attached is a patch which shows one of my failed attempts to try something to get a real/imaginary pair out of the formulas I derived. It's completely wrong, but might illustrate our problems.
Ciao
Frank Barknecht _ ______footils.org__
There's two problems with your patch. There's a sign error in one of the expr which calculates the imaginary value. The second problem (in the same expr) was a * vs / error ( you put "/ 2*$f1" instead of / 2 / $f1). Does that make a difference?
Anyway, I jiggled with it enough to get it right. Also, I put two other expr that check to make sure the poles / zeros are in fact complex rather than real.
Wait a minute.... how about just making some long expr that will sort it out correctly and supply the zeros and poles for the real case as well. I'll add that too.
Chuck p.s. I had some trouble following your expr syntax. It makes me wonder if there's something I don't know about expr?
Hallo Chuck,
thank you soo much! All is clear now. It seems I already was close or at least on the right path. The check for real roots is very nice as well, I had a bad feeling ignoring them completely and just crossing fingers in my patch.
Now it isn't difficult anymore to replace the message calculations with signal objects and make a biquad with all signal inlets, for example for Damian's dinosaurs abstractions, that are are using the Biquad Cookbook [1] formulas for various filter types.
I wonder if dealing with the cookbook includes unnecessary conversions? The cookbook calculates parameters for the transfer function (ff1, ff2, ...) from frequency, gain and Q settings, which can be used in [biquad~] directly, but must be converted to poles and zeros for elementary filters. Calculating the poles and zeros directly from frequency, gain and Q should be much simpler, wouldn't it? An Elementary-Filter-Cookbook would be a nice project.
[1] http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt
Thanks again.
Frank
Charles Henry hat gesagt: // Charles Henry wrote:
There's two problems with your patch. There's a sign error in one of the expr which calculates the imaginary value. The second problem (in the same expr) was a * vs / error ( you put "/ 2*$f1" instead of / 2 / $f1). Does that make a difference?
Anyway, I jiggled with it enough to get it right. Also, I put two other expr that check to make sure the poles / zeros are in fact complex rather than real.
Wait a minute.... how about just making some long expr that will sort it out correctly and supply the zeros and poles for the real case as well. I'll add that too.
Chuck p.s. I had some trouble following your expr syntax. It makes me wonder if there's something I don't know about expr?
Charles Henry wrote:
There's two problems with your patch. There's a sign error in one of the expr which calculates the imaginary value. The second problem (in the same expr) was a * vs / error ( you put "/ 2*$f1" instead of / 2 / $f1). Does that make a difference?
Anyway, I jiggled with it enough to get it right. Also, I put two other expr that check to make sure the poles / zeros are in fact complex rather than real.
cool, that's very nice.
however, there's some weird gain thing going on. check out the attached: notice how the left channel output (biquad~ output) drops gain, while the right channel output (handmade) does not.
any ideas how to fix this?
damian stewart | skype: damiansnz | damian@frey.co.nz frey | live art with machines | http://www.frey.co.nz
On Sat, Sep 20, 2008 at 7:09 AM, Damian Stewart damian@frey.co.nz wrote:
however, there's some weird gain thing going on. check out the attached: notice how the left channel output (biquad~ output) drops gain, while the right channel output (handmade) does not.
any ideas how to fix this?
ah! got it! Remember how I said that ff1 determines the gain? All you need to do is add a *~ block in series with the [czero~]'s and [cpole~]'s to multiply by ff1
Chuck
On Sat, Sep 20, 2008 at 12:53 PM, Charles Henry czhenry@gmail.com wrote:
On Sat, Sep 20, 2008 at 7:09 AM, Damian Stewart damian@frey.co.nz wrote:
however, there's some weird gain thing going on. check out the attached: notice how the left channel output (biquad~ output) drops gain, while the right channel output (handmade) does not.
any ideas how to fix this?
ah! got it! Remember how I said that ff1 determines the gain? All you need to do is add a *~ block in series with the [czero~]'s and [cpole~]'s to multiply by ff1
Chuck
As an afterthought, this approach does not work when ff1 is zero :) duh
But then, I don't think that it's ever useful to have ff1 be zero, is it? Also, I was wrong about the zeros, when ff1 is zero.
"Now, there's one free variable remaining, ff1. ff1 could be zero, which would give us a zero at z=0 and another real valued zero at -ff3/ff2."
There is no zero at z=0. I'm not sure about this one--but it seems as though it's impossible to have a zero at z=0?
Charles Henry wrote:
There is no zero at z=0. I'm not sure about this one--but it seems as though it's impossible to have a zero at z=0?
[shrug] i don't even know what a 'zero' is. and that's after reading several different webpages that attempted to explain biquad filtering to me. i just don't understand it.
Damian Stewart wrote:
Charles Henry wrote:
There is no zero at z=0. I'm not sure about this one--but it seems as though it's impossible to have a zero at z=0?
a zero at z=0 is a delay of 1 sample, so in:
y = a x(0) + b x(-1) + c x(-2) + d y(-1) + e y(-2)
there is a zero at z=0 if a == 0, and two zeros at z=0 if a == 0 and b == 0
(iirc)
[shrug] i don't even know what a 'zero' is. and that's after reading several different webpages that attempted to explain biquad filtering to me. i just don't understand it.
if you have a transfer function like:
(z-a)(z-b)
H(z) = g ---------- (z-c)(z-d)
then 'a' and 'b' are zeros and 'c' and 'd' are poles.
The gain+phase response at a given frequency f is given by:
H( exp(i w) ) = H( cos(w) + i sin (w) )
where w = 2 * pi * f / SR
Intuitively the closer the point exp(iw) is to a zero, the less the gain (and if it's equal to a zero the gain is zero), and the closer the point is to a pole, the greater the gain (and if it's equal to a pole, the gain is infinite - so keep pole radius strictly less than 1 !).
Claude
On Sun, Sep 21, 2008 at 8:56 AM, Claude Heiland-Allen claudiusmaximus@goto10.org wrote:
Damian Stewart wrote:
Charles Henry wrote:
There is no zero at z=0. I'm not sure about this one--but it seems as though it's impossible to have a zero at z=0?
a zero at z=0 is a delay of 1 sample, so in:
y = a x(0) + b x(-1) + c x(-2) + d y(-1) + e y(-2)
there is a zero at z=0 if a == 0, and two zeros at z=0 if a == 0 and b == 0
(iirc)
z^-1 is the unit delay operator.... so, Y(z)=z*X(z) means y(n) = x(n+1)
while this is not impossible... it's non-causal. For real-time filtering, you can't already know the sample that comes next. I've been thinking about it for a little while now. I hope we can figure it out
Chuck
Charles Henry wrote:
z^-1 is the unit delay operator.... so, Y(z)=z*X(z) means y(n) = x(n+1)
while this is not impossible... it's non-causal. For real-time filtering, you can't already know the sample that comes next. I've been thinking about it for a little while now. I hope we can figure it out
oops, i got it backwards, thanks for the correction - poles at 0 are unit delay, zeros at 0 are unit advance
but that means that rzero~ and czero~ as implemented in Pd have an extra pole at 0, so implementing a biquad~ with czero~ and cpole~ will give a 2 sample delay compared to the "real" biquad~.
snippet of source of sigrzero_perform():
{for (i = 0; i < n; i++){ t_sample next = *in1++; t_sample coef = *in2++; *out++ = next - coef * last; last = next; }}
when coef = 0, then out = in : no delay or advance, the extra pole at 0 cancels out the zero at 0
Chuck
Claude
A window or block is simply a buffer. You still only know the samples that have already passed. For most of the samples in the block, you know what sample comes next, only because it also has already passed. You don't know the sample that comes next past the boundary of the block.
Now, for example, you could program a filtering algorithm on a (sufficiently fast) FPGA with zero (not quite zero, but less than one sample) latency. That's where the causal/non-causal distinction is most clear.
On Sun, Sep 21, 2008 at 10:42 AM, Steffen Juul stffn@dibidut.dk wrote:
On 21/09/2008, at 16.02, Charles Henry wrote:
For real-time filtering, you can't already know the sample that comes next.
Isn't that what windows (=(?) blocks in Pd-lingo) is for?
Damian Stewart wrote:
[shrug] i don't even know what a 'zero' is.
a (probably a little too simple) explanation:
a zero is a valley. its depth is 0 on its deepest point. a pole is a very (very,very,...) high mountain. both things exist on the complex plane (just a plain, where all points have 2 coordinates, like any map)
you can design filters by placing (defining a pair of coordinates) poles and zeros within the unit circle (a circle with radius==1 around the origin (center) of the complex plain). the frequencys from zero to samplerate/2 are wrapped around the upper half of the unit circle. when you cut along this line through your 'landscape', the profile is the same as the graph on your filter/eq plugin.
for example the peak-filter from the pd docs is just a 'mountain' and a 'valley' which are moved on circular pathes with a certain distance to the unit circle. the closer the 'mountain' is to the unit circle, the higher the peak (on the circular profile-cutout). the frequencys coressponding to this spot are boosted. the position on the path (or the angle) defines the frequency.
so it is possible to make filters with the cpole/zero objects, without knowing too much (or anything at all...) of the math involved.
the only problem i still have, is the gain-factor. while this is easy to overcome with 'static' filters (just scale the filters input or output until it is loud enough/stops clipping), i have not yet found a method (that works without thinking (and more important: computing) too much) for time-variable-filters. ok, a table would do the trick, but thats not elegant at all...
bis denn! martin
On Mon, Sep 15, 2008 at 9:39 AM, Damian Stewart damian@frey.co.nz wrote:
Frank Barknecht wrote:
Hallo, Frank Barknecht hat gesagt: // Frank Barknecht wrote:
1*z^2 - fb1*z^1 - fb2 = 0
zb = (-fb1 +- sqrt(fb1*fb1 - 4*fb2)) / 2
Sorry, the one above is wrong, of course. Given the minus-signs in the original transfer formula it should be:
zb = (fb1 +- sqrt(fb1*fb1 + 4*fb2)) / 2
if I'm right.
see, Chuck, if Frank gets it wrong what chance do i have? ;-)
We all have problems with math ;)
This all looks good so far. zf = (-ff2 +- sqrt(ff2*ff2 - 4*ff1*ff3)) / 2*ff1 (zeros)
zb = (fb1 +- sqrt(fb1*fb1 + 4*fb2)) / 2 (poles)
I am afpd at the moment, so I can't check the documentation for cpole~/czero~ and throw in my 2 cents on the implemenation. But I can finish out the problem of finding ff1, ff2, ff3, fb1, and fb2 in terms of zf+ zf- and zb+ zb-
first, (the signs are clunky when writing them out in arthmetic, so I'll just rename them) and I'm going to write it out in a slightly different form. zf1 = zf+ = -ff2/(2*ff1) + sqrt(ff2^2/(4*ff1^2) - ff3/ff1) zf2 = zf- = -ff2/(2*ff1) - sqrt(ff2^2/(4*ff1^2) - ff3/ff1)
zb1 = zb+ = (fb1 + sqrt(fb1^2 + 4*fb2)) / 2 zb2 = zb- = (fb1 - sqrt(fb1^2 + 4*fb2)) / 2
Then, zf1 + zf2 = -ff2/ff1 (real valued) zf1 - zf2 = 2*sqrt( ff2^2/(4*ff1^2) - ff3/ff1) (could be real or imaginary)
zb1 + zb2 = fb1 (real valued) zb1 - zb2 = sqrt(fb1^2 + 4*fb2) (could be real or imaginary)
The second one is easy to solve, and we will come back to the 1st one after that.
fb1 = zb1 + zb2 fb2 = ( (zb1 - zb2)^2 - (zb1 + zb2)^2 )/4
the other has 3 variables and 2 equations. so, we make our variables into 2 groups ff2/ff1 = - (zf1 + zf2) ff3/ff1 = ( (zf1 + zf2)^2 - (zf1 - zf2)^2 )/4
Now, there's one free variable remaining, ff1. ff1 could be zero, which would give us a zero at z=0 and another real valued zero at -ff3/ff2.
Otherwise, we take ff1 to be a scaling factor, which needs to be chosen by some other condition (gain at some chosen frequency, most likely).
Chuck
-
Charles Henry -
Claude Heiland-Allen -
Damian Stewart -
Frank Barknecht -
martin brinkmann -
Steffen Juul