[PD] efficient approximation of trig functions for hi pass formula (was: could vanilla borrow iemlib's hi pass filter recipe?)

katja katjavetter at gmail.com
Wed Oct 19 19:07:09 CEST 2016

On Wed, Oct 19, 2016 at 4:25 PM, Jonathan Wilkes <jancsika at yahoo.com> wrote:
> When implemented in C, which approach takes the least amount of time
> to read, reason about, and fully comprehend?

That is an important question. Pd code is full of clever tricks and
bit hacks for dsp efficiency. What is the origin of q8_rsqrt(), why
and how does it work? What about PD_BIGORSMALL() in m_pd.h? And the
mysterious UNITBIT32 number in d_osc.c? Ideally such code should be
commented not only to denote its function (if necessary) but also to
reference the origin so you may be able to find info.

An approximation for a trig function should go in an (inline)
function, with a comment if the name can't clarify the function
sufficiently. But to fully comprehend is a different matter. Dsp code
in general takes substantial background to understand. You could
wonder why and how the approximation works, but the same question goes
for the function that it replaces.


> -Jonathan
> ________________________________
> From: katja <katjavetter at gmail.com>
> To: "pd-list at lists.iem.at" <pd-list at lists.iem.at>
> Sent: Wednesday, October 19, 2016 9:06 AM
> Subject: [PD] efficient approximation of trig functions for hi pass formula
> (was: could vanilla borrow iemlib's hi pass filter recipe?)
> Changing the thread title to reflect the new approach. Extract of the
> original thread;
> - I suggested using iemlib's hi pass filter recipe to improve
> frequency response of [hip~]
> - Christof Ressi pointed to formula in
> http://www.arpchord.com/pdf/coeffs_first_order_filters_0p1.pdf
> - this formula calculates feedback coefficient k = (1 - sin(a)) /
> cos(a) where a = 2 * pi * fc / SR
> - the filter implementation is y[n] = (x[n] - x[n-1]) * (1 + k) / 2
> +  k * y[n-1]
> - following convention in d_filter.c (and pd tilde classes in
> general), trig functions should best be approximated
> - Cyrille provided libre office linear regression result for
> (1-sin(x))/cos(x)
> Thanks for the useful infos and discussion. My 'math coach' suggested
> using odd powers of -(x-pi/2) in an approximation polynomial for
> (1-sin(x))/cos(x). The best accuracy/performance balance I could get
> is with this 5th degree polynomial:
> (-(x-pi/2))*0.4908 - (x-pi/2)^3*0.04575 - (x-pi/2)^5*0.00541
> Using this approximation in the filter formula, response at cutoff
> frequency is -3 dB with +/-0.06 dB accuracy in the required range 0 <
> x < pi. It can be efficiently implemented in C, analogous to an
> approximation Miller uses in [bp~]. So that is what I'll try next.
> Attached patch hip~-models.pd illustrates and compares filter recipes
> using vanilla objects:
> - current implementation, most efficient, accuracy +/- 3 dB
> - implementation with trig functions, least efficient, accuracy +/- 0.01 dB
> - implementation with approximation for trig functions, efficient,
> accuracy +/- 0.06 dB
> A note on efficiency: coefficients in [hip~] are only recalculated
> when cutoff frequency is changed. How important is performance for a
> function rarely called? I'm much aware of the motto 'never optimize
> early', yet I spent much time on finding a fast approximation, for
> several reasons: it's a nice math challenge, instructive for cases
> where performance matters more, and I want to respect Miller's code
> efficiency when proposing a change. Today pd is even deployed on
> embedded devices so the frugal coding approach is still relevant.
> After 20 years.
> Katja
> On Tue, Oct 18, 2016 at 10:28 AM, cyrille henry <ch at chnry.net> wrote:
>> Le 18/10/2016 à 00:47, katja a écrit :
>>> The filter recipe that Christof pointed to was easy to plug into the C
>>> code of [hip~] and works perfectly. But when looking further in
>>> d_filter.c I came across an approximation function 'sigbp_qcos()' used
>>> in the bandpass filter. It made me realize once more how passionate
>>> Miller is about efficiency. I'm not going to make a fool of myself by
>>> submitting a 'fix' using two trig functions to calculate a filter
>>> coefficient when a simple approximation could do the job. So that is
>>> what I'm now looking into, with help of a math friend: an efficient
>>> polynomial approximation for (1-sin(x))/cos(x).
>> according to libre office linear regression, for x between 0 and Pi,
>> (1-sin(x))/cos(x) is about :
>> -0.057255x³ + 0.27018x² - 0.9157x + 0.99344
>> the calc is in attachment, if you want to tune the input source or
>> precision.
>> cheers
>> c
> _______________________________________________
> Pd-list at lists.iem.at mailing list
> UNSUBSCRIBE and account-management ->
> https://lists.puredata.info/listinfo/pd-list

More information about the Pd-list mailing list