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

[Christof Ressi]

> that's where some setup code before the actual loop can become a bottleneck.

Hmmm, I forgot that in the case of [hip~] the coefficient is not calculated in the perform routine at all, so maybe efficiency is really not that important. I guess, you can't change the cutoff frequency so often that the CPU cost will actually matter. But maybe I'm wrong... 

> Gesendet: Mittwoch, 19. Oktober 2016 um 23:22 Uhr
> Von: "Christof Ressi" <christof.ressi at gmx.at>
> An: katja <katjavetter at gmail.com>
> Cc: "pd-list at lists.iem.at" <pd-list at lists.iem.at>
> Betreff: Re: [PD] efficient approximation of trig functions for hi pass formula (was: could vanilla borrow iemlib's hi pass filter recipe?)
>
> > And the mysterious UNITBIT32 number in d_osc.c?
> 
> Funny you mentioned that! I just wanted to add that instead of approximating trig functions through polynomials one could also read Pd's cosine table ( float * cos_table ) with linear interpolation like it's done in [cos~] and [osc~]. Miller does exactly this in the code for [vcf~] which - I think - is therefore part of d_osc.c and not d_filter.c.
> I'm wondering which of these methods is more efficient and/or more precise...
> 
> > A note on efficiency: coefficients in [hip~] are only recalculated
> > when cutoff frequency is changed. How important is performance for a
> > function rarely called?
> 
> It can be relevant for small blocksizes. Right now I'm working on a project where all audio is done with blocksize 1 + upsampling and that's where some setup code before the actual loop can become a bottleneck. At [block~ 1 0 8] one can get really paranoid of efficiency :-)
> 
> Christof
> 
> > Gesendet: Mittwoch, 19. Oktober 2016 um 19:07 Uhr
> > Von: katja <katjavetter at gmail.com>
> > An: "Jonathan Wilkes" <jancsika at yahoo.com>
> > Cc: "pd-list at lists.iem.at" <pd-list at lists.iem.at>
> > Betreff: Re: [PD] efficient approximation of trig functions for hi pass formula (was: could vanilla borrow iemlib's hi pass filter recipe?)
> >
> > 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.
> > 
> > Katja
> > 
> > >
> > > -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
> > >
> > >
> > 
> > _______________________________________________
> > Pd-list at lists.iem.at mailing list
> > UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
> >
> 
> _______________________________________________
> Pd-list at lists.iem.at mailing list
> UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
>