Author Topic: Calculating filter coefficient from a transfer function - literature not clear  (Read 2415 times)

0 Members and 1 Guest are viewing this topic.

Offline YansiTopic starter

  • Super Contributor
  • ***
  • Posts: 3930
  • Country: 00
  • STM32, STM8, AVR, 8051
Hello,

just got into reading a book from Udo Zölzer, DAFX Digital audio effects (link). I've been researching about programming a signal limiter/compressor and found this awesome book.

So straight into the Chapter 5, on page 98, in Figure 5.6 there are peak and RMS detector systems:



 On the bottom of that page is an equation, stating all coefficient there (AT, RT and TAV) can be calculated using this very same formula:



Details (but not enough) about that formula are provided on the next page, 99.  I got quite curious about the constant "2.2" in the exponent.  There is a reference to a [McN84] (jeez I hate this style of references!) literature, this one: BBC Digital Audio Dynamic Range Control report (link).

Relevant stuff starts at page 10, derivation of the equation is presumably on the page 11:



However, I do not see any formula with either Log or Exp within, only this one seem to be related:



It calculates the time constant of the filter: tau = 1 / 2*Pi*fc (presumably??)
If I plug in some numbers in the equation from Zölzer,  for example 10ms @48kHz fs, I get the coefficient to be A = 0.004572845.

Now if I plug this A constant into the equation from the BBC paper, I get a completely different time constant back, something like 4.55ms.

Yes, the two equations does not seem to be related at all, I was only silly trying that. But how the hell did Zölzer come up with his exponential equation, based on the reference paper?


//Edit!:   Supposed the low pass configuration in the system above is indeed equivalent to the classic single-pole filter, then the coefficient shall be
A = 1 - Exp(-2Pi*fc/fs) = 1 - Exp(-2Pi*Ts/t), shouldn't it?

The PDF seems to have gone through OCR, isn't this just a misprint, or misOCRed respectively?  :o
« Last Edit: November 30, 2019, 04:59:29 pm by Yansi »
 
The following users thanked this post: RoGeorge

Offline bson

  • Supporter
  • ****
  • Posts: 2497
  • Country: us
What is tau?  In many filters it's a parameter that controls decay, spread, number of lobes etc, which is the number of samples, the time delay, or some other aspect (e.g. # of lobes in a sinc filter), before the output is 0. If so, 2.2 might be a constant that is "a little more than nyquist".
 

Offline Wimberleytech

  • Super Contributor
  • ***
  • Posts: 1134
  • Country: us
I am not certain I can follow what your question is.  However, I did duplicate the analysis from the paper and discovered there is an error in one of their steps.  The final result for omega (-3dB) is correct, however.
 

Offline YansiTopic starter

  • Super Contributor
  • ***
  • Posts: 3930
  • Country: 00
  • STM32, STM8, AVR, 8051
Sorry if my question was not clear enough.

I think, that the equation for calculating the AT, TAV and RT coefficients is wrong. And after bit of fiddling in Octave, I can prove it is indeed wrong.

We are talking about this eqn from the Zölzer, page 98:



this eqn should give us the constant to this kind of one-pole low pass filter structure, as seen from the figures 5.6, 5.7 and others (a section of Fig 5.7):



The time constant, is the time constant of the filter, i.e as in the BBC paper:  omega = 1/tau.   Tau is either the AT or RT  time constant (attack/release time constant) of the filter element.

By plugging some numbers in Octave, I have concluded, that the equation K = 1 - Exp(-2.2 * T/t) is flawed. The measured time for the exponential to reach 63% was very incorrect.

I am certain, that this equation for a one-pole filter is correct:  K = 1 -  Exp(-2*Pi * fc/fs)
By substituing both fs (sampling frequency) with Ts (sampling period) and fc (cutoff frequency) with the time constant (which relates to fc = 1/(2*Pi*tau), we then get:

K = 1 -  Exp(-T/t)

Verifying in Octave shows a correct result.  I have put into Octave the whole structure from Fig 5.7 (Zölzer pg 98). The Octave .m file is attached, if you want to look.

I have absolutely no idea, where or how the 2.2 constant got there. It ain't right. Plot below with the corrected formula above, yields correct time constants for the filters. Plotted for  4ms attack and 50ms decay.



Also, I noted that the RMS detector in Fig 5.6 has an error in it, apart from the missing square root on the output: There should very likely be a minus sign at that summing node.




.m file in zip archive:
* AT_RT_filter.zip (0.75 kB - downloaded 43 times.)

//EDIT: FUCKING HELL, the forum randomly swaps the attachments. SHIT!!!
//EDIT2: COrrected.
//EDIT3: Added time constant description for the plot.
« Last Edit: November 30, 2019, 10:41:57 pm by Yansi »
 

Offline Wimberleytech

  • Super Contributor
  • ***
  • Posts: 1134
  • Country: us

By plugging some numbers in Octave, I have concluded, that the equation K = 1 - Exp(-2.2 * T/t) is flawed. The measured time for the exponential to reach 63% was very incorrect.


I wonder if the time constant is based on something other than 63% ?  Perhaps 90% or so???
 

Offline YansiTopic starter

  • Super Contributor
  • ***
  • Posts: 3930
  • Country: 00
  • STM32, STM8, AVR, 8051
I think a time constant of a first order system is pretty well defined. NEver heard any else.

https://en.wikipedia.org/wiki/Time_constant
 

Offline YansiTopic starter

  • Super Contributor
  • ***
  • Posts: 3930
  • Country: 00
  • STM32, STM8, AVR, 8051
If I calculate back using the 2.2 constant within that equation, the "time constant" then correcsponds to 1-Exp(-2.2) = 88.9%. Close to 90%. Maybe intentional?

Not sure, in any case the literature SHOULD BE SPECIFIC a non-standard definition of a time constant has been used for whatever reason.
 

Offline Wimberleytech

  • Super Contributor
  • ***
  • Posts: 1134
  • Country: us
If I calculate back using the 2.2 constant within that equation, the "time constant" then correcsponds to 1-Exp(-2.2) = 88.9%. Close to 90%. Maybe intentional?

Not sure, in any case the literature SHOULD BE SPECIFIC a non-standard definition of a time constant has been used for whatever reason.

Agree
 

Offline YansiTopic starter

  • Super Contributor
  • ***
  • Posts: 3930
  • Country: 00
  • STM32, STM8, AVR, 8051
Mmm... seems there is just confusion if anything. According to wiki and the cited source there, quote:


Quote
The length of each period is determined by the rate of change and the required change in gain. For more intuitive operation, a compressor's attack and release controls are labeled as a unit of time (often milliseconds). This is the amount of time it takes for the gain to change a set amount of dB or a set percentage towards the target gain. There is no industry standard for the exact meaning of these time parameters.

https://en.wikipedia.org/wiki/Dynamic_range_compression#Attack_and_release

 >:(
 

Offline YansiTopic starter

  • Super Contributor
  • ***
  • Posts: 3930
  • Country: 00
  • STM32, STM8, AVR, 8051
In this document from Powersoft, they state that the reaction times are from the beginning of the event till the return to the previous state  - so more like the 90-95% in Zölzer, rather than a standard definition of 63% tau.


See page 2, sections 1.2, 1.3: http://www.powersoft-audio.com/en/docman/804-how-to-setup-limiters/file

Okay, it may make sense after all for this application. Still not sure if this is the "industry accepted" behavior, but at least two sources are somewhat similar (Zölzer, Powersoft).

I tried to look for this also in some user manuals from Alesis and DBX compressors, no luck. Manuals for too much of a n00b user to be this advanced and comprehensive, as the Powersoft manual is.
« Last Edit: December 02, 2019, 07:12:13 pm by Yansi »
 

Offline RoGeorge

  • Super Contributor
  • ***
  • Posts: 7012
  • Country: ro
Thanks for the book!
Just read a little the IIR/FIR chapter and loved it.  I like the stile it is written.   :-+

Offline YansiTopic starter

  • Super Contributor
  • ***
  • Posts: 3930
  • Country: 00
  • STM32, STM8, AVR, 8051
So, after playing a bit and learning to script in Octave, I have finally managed to make the attack/release aware envelope detection algorithm to work.

It is not that difficult, after all. Bud mind you, I am just a DSP (and Octave) beginner, so I am pretty satisfied with the result  ;D

Code: [Select]
function pk = peakAtRt(in, AT, RT)
  persistent lastpk = 0;
  a = abs(in) - lastpk;
  if (a < 0) a = 0; endif;
  pk = AT*a + lastpk*(1-RT);
  lastpk = pk;
endfunction;



 

Offline YansiTopic starter

  • Super Contributor
  • ***
  • Posts: 3930
  • Country: 00
  • STM32, STM8, AVR, 8051
...and do even some rudimentary compression  ;D

 

Offline rhb

  • Super Contributor
  • ***
  • Posts: 3516
  • Country: us
For an appropriate range of values, calculate the filter element transfer functions organized such that the columns in the A matrix are the candidate source population transfer functions.  Choose the exact points in your parameter space of interest randomly. It is extremely important that the points in the parameter space of A show no correlations.

If you then solve Ax=y where y is the desired transfer function,  A is a basis pursuit dictionary  and x is a sparse  L1 solution using linear programming  you'll have the provably optimum terms. 

If you have parts with measured values such that you can select which parts are used together in assembly you can substantially improve system performance for the price of a few compute cycles.  With a 100% incoming test on reel packed parts the data are available and the only constraint is needing random pickup from the tape.

You seek the optimal combination and you are overwhelmingly likely to get an answer which is provably optimal.  But sometimes you don't get an answer.  Not often, but sometimes.

But if you can cast your problem as the product of a bunch of unknown transfer functions.  So you make up a large enough random sample of the transfer functions, poles and zeros.  Linear programming doesn't solve problems like that unless you take the logarithm in which case the transfer function is the sum of terms and a linear problem.

The arbitrary precision simplex solver in GLPK is spectacularly good at doing these.  And glpsol and GMPL are a very convenient interface.

The classical way of doing this is the Wiener filter which solves Ax=y using an L2 norm.   That actually works very well.  It has the disadvantage that L2 smears the answer and produce unstable filters. L1 won't do that if you restrict membership in the dictionary to stable poles and zeros.  But L2 requires much less compute effort.   So it's by no means obsolete.  Just not the only game in town as when L2 ruled because of CPU and memory cost.

Have Fun!
Reg
 

Offline YansiTopic starter

  • Super Contributor
  • ***
  • Posts: 3930
  • Country: 00
  • STM32, STM8, AVR, 8051
Eehhh.... you sure the reply wasn't meant to a different thread?  :-//
 

Offline rhb

  • Super Contributor
  • ***
  • Posts: 3516
  • Country: us
You are concerned with providing the best match to an arbitrary transfer function  under certain constraints as to performance and cost aren't you?

Matched filter design is routine in seismic and most other DSP environments.  What would you like the spectrum to look like?  Do you just want me to search for that recordings containing that voice spectrum?

Wiener developed all of this around 1940 and later published as "The Smoothing, Interpolation and Extrapolation of Stationary Time Series." in 1949.  In it he presents the convolutional model, one of  the most basic models in mathematical physics.

Candes and Donoho are as big a step change as Wiener and Shannon were.

If you want to look like a genius, learn to formulate and solve problems such as the one you posed via linear programming.  You and all around you will be amazed.  It is the wildest stuff I've ever seen. 

Have Fun!
Reg
 

Offline YansiTopic starter

  • Super Contributor
  • ***
  • Posts: 3930
  • Country: 00
  • STM32, STM8, AVR, 8051
Well, I am not trying to best match a transfer function, hence why I don't understand you. I just try understanding, how to make a signal limiter in a DSP (mainly for speaker rms power protection and to prevent peaks from saturating the power amplifier)

And after toying with it for a while, I am finding more and more issues with the systems provided in Zölzer. He has likely just copied from the BBC paper, without even trying to verify.  :-\

I am pretty sure now the AT/RT peak detector system is not working correctly. See below. I am talking about this system, which is supposed to provide signal peak value detection, while respecting a set attack and release time:



More thorough testing with simple waveforms, where the behavior shall be very predictable revealed, that the peak value given in steady state is incorrect:



After poking at it a bit, I see that by providing very short attack times (AT constant is close to 1), the result is correct. So it is also when RT constant is close to 0.
Now thinking about how to patch the system, I am thinking about the feedback path to the |x(n)| subtracting node shall be scaled by some amount.
« Last Edit: December 03, 2019, 10:45:35 pm by Yansi »
 

Offline RoGeorge

  • Super Contributor
  • ***
  • Posts: 7012
  • Country: ro


I got quite curious about the constant "2.2" in the exponent.

Digital came quite late in audio processing.  Before that, it was all analog, so without following too close what is trying to implement there, for sure it's some former analog behavior emulated digitally.  That given said, anything with a time constant in the analog world back then was done with RC circuits.  Looking at the time constant of an RC circuit:

https://en.wikipedia.org/wiki/RC_time_constant
Quote from: Wikipedia
time (10% to 90%) tr ≈ 2.2 τ ≈ 0.35/fc

My best guess is that's where from it was coming the 2.2 magic number you were asking about.
« Last Edit: December 04, 2019, 01:37:40 am by RoGeorge »
 
The following users thanked this post: Wimberleytech


Share me

Digg  Facebook  SlashDot  Delicious  Technorati  Twitter  Google  Yahoo
Smf