[PD] OT : using libre office for data regression; WAS Re: efficient approximation of trig functions for hi pass formula (was: could vanilla borrow iemlib's hi pass filter recipe?)

katja katjavetter at gmail.com
Thu Oct 20 01:07:59 CEST 2016


On Wed, Oct 19, 2016 at 6:32 PM, cyrille henry <ch at chnry.net> wrote:
>
>
> Le 19/10/2016 à 17:34, katja a écrit :
>>
>> Thanks for your update Cyrille. This seems a useful approach to find
>> approximations. I should really learn how to do that with libre
>> office. Do you have tutorial links? So far I found this:
>>
>> https://newtonexcelbach.wordpress.com/2015/06/28/using-linest-for-non-linear-curve-fitting-examples-hints-and-warnings/.
>
>
> I  learn to do this yesterday...
> Most informations I found on the web was outdated, but it's very simple.
> for short :
> create 2 column with X/Y data.
> select the data.
> click insert / diagram. select X/Y dispersion, in order to draw a X/Y graph
> of your data.
> then double click on the diagram and click on the curve to select it.
> then right click and "insert a tendency curve"
> then enjoy all option available.
> (there is no "odd only polynomial", but even value was almost 0 in the last
> example).
>
> Don't forget to click "draw the equation"
> The only problem is that i did not find a way to cut/paste the equation...
>
> (my libre office is in french, so menu translation may not be accurate)
>
>
>> For the filter recipe the curve must go through coordinates [0 1] and
>> [pi -1], to get the expected behavior when cutoff frequency is set to
>> DC or Nyquist. [hip~ 0] should not block DC, which is only possible
>> when the coefficient is exactly 1 like you get it with
>> (1-sin(0))/cos(0). I tuned my factors 'by hand' to do so. That may be
>> possible with your libre office result as well.
>
> yes, in my example the result is slightly different from 1 / -1.
> i tried again, adding about 20 points for X=0 and X=1 to get more weight for
> this coordinate.
> the coef where better, but not perfect.
>
> still it's a good starting point for manual adjustment.

That's what I thought, avoiding the initial guess iterations saves a
lot of time. Can you post your adjusted parameters? Being an
approximation the function will never be perfect. We can select the
best frequency response compromise.

>
>
> However, when done in
>>
>> double precision the calculation will give slightly different result.
>> Above or below 1, I don't know yet. In any case this has to be
>> considered when writing the C.
>
> shouldn’t it be clipped? (what is a filter with negative frequency, or
> frequency above Nyquist?)

Yes coefficients below -1 and above 1 must be clipped. But what if it
doesn't reach 1 at DC or -1 at Nyquist in double precision because of
the rounding difference. The behavior at those extremes would be
incorrect. Not that it matters now, we don't have double precision pd.
But code written today should be double-precision-proof, you never
know what will happen tomorrow.

>
> cheers
> c
>
>>
>> Katja
>>
>> On Wed, Oct 19, 2016 at 4:18 PM, cyrille henry <ch at chnry.net> wrote:
>>>
>>> since you manually adjust the coefficient, i wanted to see the difference
>>> with coefficient adjusted/optimised by a computer.
>>> I update the calc  using X from -0.5 to 0.5 with X = (x-0.5)*Pi
>>> The coefficient i get are really different from yours (when considering
>>> the
>>> Pi factor between them).
>>>
>>> I update you patch with this coefs.
>>> I don't understand why, but your coefs works better for low frequency.
>>> so, nothing better here, but I still think it's worth sharing, in case
>>> i'm
>>> not the only one wondering...
>>>
>>>
>>> also, there was a typo in my original mail, it's not a linear regression,
>>> but a polynomial regression (since the result was obviously polynomial)
>>>
>>> cheers
>>> c
>>>
>>>
>>> Le 19/10/2016 à 15:06, katja a écrit :
>>>>
>>>>
>>>> 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
>>
>
> _______________________________________________
> 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