Discussion:
[music-dsp] How to calculate Leaky Integrator cutoff frequency?
Steven Cook
2009-08-22 08:40:22 UTC
Permalink
Hi,

My subject says it all really - what is the correct way (if any) of
calculating the cutoff frequency of a Leaky Integrator. I'm currently using
the following formula, but I've no idea where I got it from or if it is
correct, so I'm starting to doubt that it is...

cutoff = pow (e, -fc / sampleRate * PI * 2.0);

...where 'fc' is the desired cutoff frequency and 'e' is the well-known
mathematical constant.

Regards,

Steven Cook.
http://www.spcplugins.com/
+44(0)1271 867288
Bogac Topaktas
2009-08-22 16:25:03 UTC
Permalink
It is an approximation for the 3 dB cutoff frequency of an exponential averager (or leaky integrator). As the cutoff frequency increases, the actual 3dB cutoff point deviates proportionally (for instance Fc: 5000 Hz, fs: 44000 Hz, -3db point: 5213 Hz).

alpha = 1 - e^(-(Fc / Fs)*2*pi) // 1 - pow (e, -(Fc / Fs) * PI * 2.0);

where the difference equation is:

y(n) = alpha * x(n) + (1 - alpha) * y(n-1)

or re-arrange to save one multiplier:

y(n) = y(n-1) + alpha * ( x(n) - y(n-1) )

For more information see:

http://blog.prosig.com/2003/04/28/data-smoothing-rc-filtering-and-exponential-averaging/

http://en.wikipedia.org/wiki/Exponential_smoothing

http://www.dsprelated.com/showarticle/72.php

http://www.ibrtses.com/embedded/exponential.html

-------- Original Message --------
Sent: Saturday, August 22, 2009 11:41 AM
Subject: SPAM-LOW: [music-dsp] How to calculate Leaky Integrator cutoff frequency?
Hi,
My subject says it all really - what is the correct way (if any) of
calculating the cutoff frequency of a Leaky Integrator. I'm currently using
the following formula, but I've no idea where I got it from or if it is
correct, so I'm starting to doubt that it is...
cutoff = pow (e, -fc / sampleRate * PI * 2.0);
...where 'fc' is the desired cutoff frequency and 'e' is the well-known
mathematical constant.
Regards,
Steven Cook.
http://www.spcplugins.com/
+44(0)1271 867288
--
subscription info, FAQ, source code archive, list archive, book reviews, dsp links
http://music.columbia.edu/cmc/music-dsp
http://music.columbia.edu/mailman/listinfo/music-dsp
robert bristow-johnson
2009-08-22 21:36:04 UTC
Permalink
it's nice to see music-dsp happening again. i thought it was falling
silent.
Post by Bogac Topaktas
y(n) = alpha * x(n) + (1 - alpha) * y(n-1)
y(n) = y(n-1) + alpha * ( x(n) - y(n-1) )
in my opinion, the simplest way to look at this is to recognize that
a leaky integrator is the simpliest 1-pole digital low-pass filter
(not to be confused with the 1-pole, 1-zero LPF that you would get
from bilinear transforming a simple analog 1-pole LPF). here, the
pole, p, is 1-alpha

y[n] = (1-p)*x[n] + p*y[n-1]

Y(z) = (1-p)*X(z) + p*(z^-1)*Y(z)


Y(z)/X(z) = H(z) = (1 - p) / (1 - p*z^-1)

= (1-p)*z / (z - p)

H(e^jw) = (1-p)/(1 - p*cos(w) + j*sin(w))

the magnitude squared frequency response:

|H(e^jw)|^2 = (1-p)^2 / ( (1 - p*cos(w))^2 + (sin(w))^2 )

= (1-p)^2 / (1 + p^2 - 2*p*cos(w))

the -3dB corner frequency, w0, (a.k.a. "half-power frequency") occurs
when |H(e^jw)|^2 = 1/2 or the denominator of |H(e^jw)|^2 is twice as
large as the numerator:


1 + p^2 - 2*p*cos(w0) = 2*(1-p)^2 = 2*(1 - 2*p + p^2) = 2 - 4*p
+ 2*p^2

cos(w0) = -(1 - 4*p + p^2)/(2*p) = 2 - (1/2)*(p + 1/p) = 1 -
((1-p)^2)/(2*p)

(choose your favorite form)

then w0 is


w0 = 2*pi*f0/Fs = arccos( 2 - (1/2)*(p + 1/p) )

= arccos( 2 - cosh(log(p)) )


plug in 1-alpha in for p, and it looks slightly less pretty.


you can see that since the magnitude to the argument to the arccos()
function must be no greater than 1, then


1 <= cosh(log(p)) = (1/2)*(p + 1/p) <= 3

or (remember that log(p) < 0 so we take the left-hand arccosh)

0 <= -log(p) <= arccosh(3) = log(3 + sqrt(3^2-1)) =
1.762747174

or

e^(-arccosh(3)) = 0.171572875 <= p <= 1 .

we know that p < 1 for stability, i think the meaning of restricting
0.171572875 <= p (clearly p can be set closer to zero and you still
have a working 1-pole filter) is that there is no -3 dB frequency if
p is lower.


--

r b-j ***@audioimagination.com

"Imagination is more important than knowledge."
Didier Dambrin
2009-08-23 00:58:19 UTC
Permalink
Not that it's any important in practice, but something has always puzzled me
with the dB scale.
Too often I read "+6dB" talking about a 2x gain. But I like square things
thus I never quite liked that dB scale where +6dB means *1.995something.

Example: the definition of brown noise says -6dB/octave. But if I was to
generate brown noise, I'd probably rather divide by 2 per octave, it sounds
more "logical", as I wouldn't see where those x1.995/octave would come from.

So, the dB scale being 10^(x/20) & not 2^(x/6), is it only for historical
reasons or something I'm missing?
& is it common to interpret dB values as 2^(x/6)? Even just to display
better rounded values to users?

& in fact, why aren't people using a simple 2^x linear scale for gain, just
as for octaves?
Richard Dobson
2009-08-23 01:30:16 UTC
Permalink
Post by Didier Dambrin
Not that it's any important in practice, but something has always puzzled me
with the dB scale.
Too often I read "+6dB" talking about a 2x gain. But I like square things
thus I never quite liked that dB scale where +6dB means *1.995something.
That's just a matter of where you do your casual rounding. It is just
as likely you will divide (or multiply) by 2, and the proper dB value
(20 * log10(x)) is +-6.02. Nobody wants so say "six point oh two" all
the time, unless you want to be thought of as a COMPLETE NERD, so we say
"sixdeebee" and leave it at at that. Which is still pretty nerdy. Better
to say "a fifth" than "702 Cents". Even though that is what an
equal-tempered fifth is.
Post by Didier Dambrin
Example: the definition of brown noise says -6dB/octave.
Engineers will say 20dB/decade. Take your pick!

But if I was to
Post by Didier Dambrin
generate brown noise, I'd probably rather divide by 2 per octave, it sounds
more "logical", as I wouldn't see where those x1.995/octave would come from.
The dB range is nothing more than a convenience, with respect to our
ears. We hear logarithmically, and -60dB (colloquially speaking) is the
"threshold of silence", rather more convenient and expressive than
saying "point oh oh one". We need ~something that (a) compares things
and (b) give us a tidy scale with which to discuss the huge dynamic
range we can hear (and measure). It's a convention; like feet and
inches. When you scale by 0.05, is that a lot? - not very much ? medium?
How much more is it that 0.03?
Post by Didier Dambrin
So, the dB scale being 10^(x/20) & not 2^(x/6), is it only for historical
reasons or something I'm missing?
Been around a very long time. See no reason to change it. Though you
~could~.
Post by Didier Dambrin
& is it common to interpret dB values as 2^(x/6)? Even just to display
better rounded values to users?
& in fact, why aren't people using a simple 2^x linear scale for gain, just
as for octaves?
We can't really deal with more than about ten octaves; and a semitone is
quite a large difference sometimes. But with intensity, the range is
(hand-waving) 0.0001 to 10000 (or so), and the smallest difference we
can detect is pretty big - about a dB. Golden-ears will claim a half-dB.

The electrical engineers will no doubt have their own perspective on it.
All about comparing voltages v powers, etc.

Too late at night to be more lucid, sorry!

Richard Dobson
Didier Dambrin
2009-08-23 02:17:10 UTC
Permalink
Well you'd be surprised. For years I've been asked for "more precision" in
my volume faders or meters, people wanna see all dB decimals. I think people
like rounded numbers in general (I thought that was a programmer's
disorder), no one can stand "4.11" when it can be a rounded 4.
So maybe it would have been better to use a *6 dB scale, where +1 doubles
the amplitude.
Post by Richard Dobson
But with intensity, the range is
(hand-waving) 0.0001 to 10000 (or so), and the smallest difference we
can detect is pretty big - about a dB. Golden-ears will claim a half-dB.
The electrical engineers will no doubt have their own perspective on it.
All about comparing voltages v powers, etc.
Too late at night to be more lucid, sorry!
Richard Dobson
Richard Dobson
2009-08-23 07:25:50 UTC
Permalink
Post by Didier Dambrin
Well you'd be surprised. For years I've been asked for "more precision" in
my volume faders or meters, people wanna see all dB decimals. I think people
like rounded numbers in general (I thought that was a programmer's
disorder), no one can stand "4.11" when it can be a rounded 4.
So maybe it would have been better to use a *6 dB scale, where +1 doubles
the amplitude.
Hmm, not really - as I wrote, the smallest difference of level most
people can distinguish (somewhat frequency-dependent) is 1dB. Digital
volume control chips are usually quantized in half-dB steps (and plenty
of audio engineers will say that that is the minimum acceptable). I
don't know if my Denon DM30 has that; but I can hear each step as the
volume is changed. Of course, what the typical punter thinks is "half as
loud" is very much cruder than that - I don't know the references, but
IIRC listening tests showed a wide variation in such a subjective thing.
A half-dB step is educational.

Sounds like you are really looking for a VU meter scale - -20 to +3. OK
for overall power monitoring, but not for getting the hi-hat to sit just
right in the mix.

I lopoked at good-ol' wikipedia again, and this says it most succinctly:
"A change in power ratio by a factor of 10 is a 10 dB change. ". Most
of the time, that is what we are interested in, not voltages but powers
(ratio of powers is 10*log10(p1/p2)); so the decibel scale gets you a
direct map of what is actually going on, power-wise. There could be a
case, I suppose, for nominating a 3dB change ("half-power point") as a
unit. but mastering engineers routinely set EQ changes of fractional dB.
Funnily enough, 1.5dB (or thereabouts) seems a popular increment.

Besides, your +1 scale sounds dangerously close to the guitar amp
discussion in "This is Spinal Tap" : "This goes up to eleven!".

Richard Dobson
Didier Dambrin
2009-08-23 13:45:37 UTC
Permalink
Yes but that's what I like about it: often a gain knob must allow increase &
reduction, but generally you allow more space for reduction.
So I like the idea of a 0..11 gain, 10 being the 0dB point.
At the same time I'd never make a gain control linear in dB's since it'd
have to go down to -infinite. But I like the 0..10..11 as -inf..0dB..6dB.
I in fact often calibrate my volume controls to -inf..0dB..6dB but more in a
1..4..5 (then 1..10..12) range, unless I need to allow more gain.
I like this 75-25% split, which also roughly matches the usual "100" MIDI
7bit default for velocity or volume. (I seem to remember the top 25% was
around +5dB in the MIDI group standardization)

So gain controls in software are (I assume) never linear in power nor in dB,
yet we have no standard linear scale that the user would understand.
Post by Richard Dobson
Besides, your +1 scale sounds dangerously close to the guitar amp
discussion in "This is Spinal Tap" : "This goes up to eleven!".
Richard Dobson
Richard Dobson
2009-08-23 15:24:50 UTC
Permalink
Post by Didier Dambrin
Yes but that's what I like about it: often a gain knob must allow increase &
reduction, but generally you allow more space for reduction.
So I like the idea of a 0..11 gain, 10 being the 0dB point.
At the same time I'd never make a gain control linear in dB's since it'd
have to go down to -infinite.
What's the problem there" At the end of the day you have a finite
word-length. No need to compute anything less than -144dB! Drawing -inf
on a display is just that - a picture.


But I like the 0..10..11 as -inf..0dB..6dB.

That's (among other things) a calibration issue (and an issue of where
on that scale you put 0dBFS).

It depends on whether compatibility with industry practice is important
to you. See e.g. the Bob Katz K-metering system, which systemizes
digital headroom ~and the fader law for different requirements:

http://www.digido.com/level-practices-part-2-includes-the-k-system.html


Richard Dobson
Didier Dambrin
2009-08-23 16:36:06 UTC
Permalink
But IMHO a software gain control should always go down to 0. For the
user -144dB would probably fine, but with floating point mixing, it makes a
difference whether the signal is inaudible or really silent. If it's really
silent you can skip further processing as an optimization.
And even if you stopped at -144dB, you still wouldn't use a linear scale in
dB for a gain control, as half of the control's range wouldn't be much
useful.
Post by Richard Dobson
Post by Didier Dambrin
Yes but that's what I like about it: often a gain knob must allow increase &
reduction, but generally you allow more space for reduction.
So I like the idea of a 0..11 gain, 10 being the 0dB point.
At the same time I'd never make a gain control linear in dB's since it'd
have to go down to -infinite.
What's the problem there" At the end of the day you have a finite
word-length. No need to compute anything less than -144dB! Drawing -inf
on a display is just that - a picture.
But I like the 0..10..11 as -inf..0dB..6dB.
That's (among other things) a calibration issue (and an issue of where
on that scale you put 0dBFS).
It depends on whether compatibility with industry practice is important
to you. See e.g. the Bob Katz K-metering system, which systemizes
http://www.digido.com/level-practices-part-2-includes-the-k-system.html
Richard Dobson
Vadim Zavalishin
2009-08-24 11:00:53 UTC
Permalink
Hi Didier and Richard

Thought, I'd drop in some thoughts.
Post by Richard Dobson
Hmm, not really - as I wrote, the smallest difference of level most
people can distinguish (somewhat frequency-dependent) is 1dB. Digital
volume control chips are usually quantized in half-dB steps (and plenty
of audio engineers will say that that is the minimum acceptable). I
don't know if my Denon DM30 has that; but I can hear each step as the
volume is changed.
I would guess, that the 1dB minimum audible difference is referring to the
context of two signals being played one after another with some gap between
them. If you quickly change a level of a continuously playing signal you
might be able to recognize a smaller change. Further on, I would expect that
in the case of mixing several tracks together, smaller gain changes for a
given track might be able to produce an audible difference (not neccesarily
explicitly perceived as a loudness change!) due to psychoacoustic effects of
the signal interaction. Plus, the subtleties are very important in the mix,
and audio engineers should have a good ear for that.

Regards,
Vadim
--
Vadim Zavalishin
Senior Software Developer | R&D

Tel +49-30-611035-0
Fax +49-30-611035-2600

NATIVE INSTRUMENTS GmbH
Schlesische Str. 28
10997 Berlin, Germany
http://www.native-instruments.com

Registergericht: Amtsgericht Charlottenburg
Registernummer: HRB 72458
UST.-ID.-Nr. DE 20 374 7747
Geschaeftsfuehrung: Daniel Haver (CEO), Mate Galic
Steven Cook
2009-08-23 08:41:40 UTC
Permalink
Post by Bogac Topaktas
y(n) = alpha * x(n) + (1 - alpha) * y(n-1)
y(n) = y(n-1) + alpha * ( x(n) - y(n-1) )
Thanks for these answers. My code (for a BLIT-based synthesizer) looks more
like:

y(n) = y(n-1) * alpha + x(in)

Does that make a difference?

Steven Cook.
http://www.spcplugins.com/
Lubomir I. Ivanov
2009-08-23 16:25:50 UTC
Permalink
----- Original Message -----
Post by Steven Cook
Post by Bogac Topaktas
y(n) = alpha * x(n) + (1 - alpha) * y(n-1)
Thanks for these answers. My code (for a BLIT-based synthesizer) looks more
y(n) = y(n-1) * alpha + x(n)
Does that make a difference?
for one pole (p0) and requirement of unity gain:

p0 = 0.7 (this is a filter that has a zero at z = 0)
--
a0 = g (gain) = 0.3
b0 = 1
b1 = -p0 = -0.7

difference equation:
y[n] = a0*x[n-1] - b1*y[n-1] = 0.3 * x[n-1] + 0.7 * y[n-1]

--------------------
but there is a relation:
(a0+...+a[n]) = (b0+...+b[n])

if we add a zero (z0) to the z-plane:

p0 = 0.7
z0 = -0.5
--
coefficient calculation:
a0 = g;
a1 = -g*z0 = g*0.5
b0 = 1
b1 = -p0 = -0.7

from unity gain requirement, (a0 + a1) must be equal to (b0 + b1 = 1 -
0.7) or:
a0 + a1 = 0.3;
g + g*0.5 = 0.3
g = 0.2;

therefore:
a0 = 0.2;
a1 = 0.1;

and the difference equation is:
y[n] = a0*x[n] + a1*x[n-1] - b1*y[n-1] = 0.2*x[n] + 0.1*x[n-1] +
0.7*y[n-1];


-----------
regards
lubomir
robert bristow-johnson
2009-08-23 20:00:01 UTC
Permalink
Post by Steven Cook
Post by Bogac Topaktas
y(n) = alpha * x(n) + (1 - alpha) * y(n-1)
y(n) = y(n-1) + alpha * ( x(n) - y(n-1) )
Thanks for these answers. My code (for a BLIT-based synthesizer) looks more
y(n) = y(n-1) * alpha + x(in)
Does that make a difference?
only in nomenclature and in a constant gain that does not affect the
-3 dB frequency because what we *mean* by the "-3 dB frequency" is
the frequency where the gain is 3 dB lower than at DC, no matter what
the gain there is.

your "alpha" is precisely the same as "p" in what i wrote. but,
using your alpha, the model i describe is


y[n] = (1-alpha)*x[n] + alpha*y[n-1]

the only difference is that my input is scaled by (1-alpha) and yours
is not. my DC gain is 1 (or 0 dB). that means that your DC gain is
1/(1-alpha) which is bigger than 1. the -3 dB frequency is still the
same:


w0 = 2*pi*f0/Fs = arccos( 2 - (1/2)*(alpha + 1/alpha) )

= arccos( 2 - cosh(log(alpha)) )

if your alpha is less than e^(-arccosh(3)) = 0.171572875, then you
don't have a -3 dB frequency (the LPF attenuation never exceeds 3
dB). if your alpha is equal to 1, you have a non-leaky integrator
and if your alpha is greater than 1, your LPF will blow up.

one last notational issue;\: i would recommend sticking with the now
common notation that you see in DSP and discrete-filtering texts and
most of the IEEE (or AES) lit. for discrete-time or discrete-
frequency sequences, we use brackets:

x[n], y[n], h[n], X[k], Y[k], H[k]

for continuous-time and continuous-frequency functions, we use parenths:

x(t), X(f), or maybe X(omega) or X(w) or sometimes
(for consistency with the Laplace transform), X(jw)

in older textbooks, they used the standard subscript

x or x_n instead of "x[n]"
n

but sometimes, in continuous-time signals, we would get a *vector* of
signals

{ x_1(t), x_2(t), x_3(t), ... x_m(t) }

where, if this was converted to discrete time, it was cumbersome to
have two subscripts; one for which element in the vector and another
for discrete time. so the present notation was evolved

{ x_1[n], x_2[n], x_3[n], ... x_m[n] }

and widely adopted. of course this notation doesn't work with
MATLAB, but they screw up the x[0] issue anyway (the bastards).

the main reason is that if you see an argument in brackets or as a
subscript, we're talking about a discrete sequence where the argument
better the hell be an integer rather than a continuous function where
the argument is not restricted to or generally assumed to be an integer.

--

r b-j ***@audioimagination.com

"Imagination is more important than knowledge."
Bogac Topaktas
2009-08-22 19:53:42 UTC
Permalink
Following the derivations in link 1 and applying a quick and dirty fix, I
have arrived at a formula which vastly minimizes the deviation at higher
cutoff frequencies:

h = 1 / Fs and T = 1 / 2 * pi * Fc therefore:

T / T + h = fs / ( 2 * pi * Fc + Fs ) = 1 - alpha

This gives almost equal amount of deviation in the opposite way so combine
them together(with averaging) to minimize the deviation:

temp = 2.0 * pi * fc;
alpha1 = e^( -( temp / fs ) )
alpha2 = fs / ( temp + fs )
alpha = 1.0 - ( ( alpha1 + alpha2 ) * 0.5 )

again, where the difference equation is:

y(n) = alpha * x(n) + (1 - alpha) * y(n-1)

or re-arrange to save one multiplier:

y(n) = y(n-1) + alpha * ( x(n) - y(n-1) )

Use at your own risk!

P.S. This is not the proper way of performing applied mathematics! It is
just a side way(*) that you have to take while trying to quickly solve a
problem just before going out on Saturday Night:)

(*)
http://ocw.mit.edu/OcwWeb/Mathematics/18-098January--IAP--2008/CourseHome/

-------- Original Message --------
Sent: Saturday, August 22, 2009 7:26 PM
To: "A discussion list for music-related DSP"
Subject: SPAM-LOW: Re: [music-dsp] How to calculate Leaky Integrator
cutoff frequency?
It is an approximation for the 3 dB cutoff frequency of an exponential
averager (or leaky integrator). As the cutoff frequency increases, the
actual 3dB cutoff point deviates proportionally (for instance Fc: 5000 Hz,
fs: 44000 Hz, -3db point: 5213 Hz).
alpha = 1 - e^(-(Fc / Fs)*2*pi) // 1 - pow (e, -(Fc / Fs) * PI * 2.0);
y(n) = alpha * x(n) + (1 - alpha) * y(n-1)
y(n) = y(n-1) + alpha * ( x(n) - y(n-1) )
http://blog.prosig.com/2003/04/28/data-smoothing-rc-filtering-and-exponentia
l-averaging/
http://en.wikipedia.org/wiki/Exponential_smoothing
http://www.dsprelated.com/showarticle/72.php
http://www.ibrtses.com/embedded/exponential.html
-------- Original Message --------
Sent: Saturday, August 22, 2009 11:41 AM
Subject: SPAM-LOW: [music-dsp] How to calculate Leaky Integrator
cutoff frequency?
Hi,
My subject says it all really - what is the correct way (if any) of
calculating the cutoff frequency of a Leaky Integrator. I'm currently
using
the following formula, but I've no idea where I got it from or if it is
correct, so I'm starting to doubt that it is...
cutoff = pow (e, -fc / sampleRate * PI * 2.0);
...where 'fc' is the desired cutoff frequency and 'e' is the well-known
mathematical constant.
Regards,
Steven Cook.
http://www.spcplugins.com/
+44(0)1271 867288
--
subscription info, FAQ, source code archive, list archive, book
reviews, dsp links
http://music.columbia.edu/cmc/music-dsp
http://music.columbia.edu/mailman/listinfo/music-dsp
--
subscription info, FAQ, source code archive, list archive, book reviews,
dsp links
http://music.columbia.edu/cmc/music-dsp
http://music.columbia.edu/mailman/listinfo/music-dsp
Bogac Topaktas
2009-08-23 14:01:49 UTC
Permalink
Post by Steven Cook
Thanks for these answers. My code (for a BLIT-based synthesizer) looks
more
Post by Steven Cook
y(n) = y(n-1) * alpha + x(in)
Does that make a difference?
The gain at DC (i.e. 0 Hz) would not be 0 dB but much higher (inversely
proportional to corner frequency) depending on the value of alpha
coefficient.

-------- Original Message --------
Post by Steven Cook
Sent: Sunday, August 23, 2009 11:42 AM
To: "A discussion list for music-related DSP"
Subject: SPAM-LOW: Re: [music-dsp] How to calculate Leaky Integrator
cutoff frequency?
Post by Steven Cook
Post by Bogac Topaktas
y(n) = alpha * x(n) + (1 - alpha) * y(n-1)
y(n) = y(n-1) + alpha * ( x(n) - y(n-1) )
Thanks for these answers. My code (for a BLIT-based synthesizer) looks
more
Post by Steven Cook
y(n) = y(n-1) * alpha + x(in)
Does that make a difference?
Steven Cook.
http://www.spcplugins.com/
--
subscription info, FAQ, source code archive, list archive, book reviews,
dsp links
Post by Steven Cook
http://music.columbia.edu/cmc/music-dsp
http://music.columbia.edu/mailman/listinfo/music-dsp
Loading...