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@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
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@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@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
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....
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. 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.
Katja
On Wed, Oct 19, 2016 at 4:18 PM, cyrille henry ch@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@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@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
Pd-list@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
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....
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.
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?)
cheers c
Katja
On Wed, Oct 19, 2016 at 4:18 PM, cyrille henry ch@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@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@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
Pd-list@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
Pd-list@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
On Wed, Oct 19, 2016 at 6:32 PM, cyrille henry ch@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....
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@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@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@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
Pd-list@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
Pd-list@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
Pd-list@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
Le 20/10/2016 à 01:07, katja a écrit : [...]
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.
libre office say :
-1.5743519661 X -1.1794365231 X^3 -2.0916936873 X^5 (for X between -0.5 / 0.5)
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.
since we know that aX + bX³ + cX⁵ should be 1 for X = -0.5, we can automatically adjust "a" to: a = 1 - b*0.5³ - c*0.5⁵ using only b and c from the polynomial regression.
(since the function is odd, it will also work for X = 0.5)
I did not try, but apart from rounding error, it should satisfy the extreme condition whatever precision is used, as long as "a" is computed in the C code.
(the pi value I used in libreoffice is not double precision accurate anyhow).
cheers c
When implemented in C, which approach takes the least amount of time to read, reason about, and fully comprehend?
-Jonathan
From: katja <katjavetter@gmail.com>
To: "pd-list@lists.iem.at" pd-list@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@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@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
On Wed, Oct 19, 2016 at 4:25 PM, Jonathan Wilkes jancsika@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@gmail.com To: "pd-list@lists.iem.at" pd-list@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@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@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
On Wed, Oct 19, 2016 at 1:07 PM, katja katjavetter@gmail.com wrote:
On Wed, Oct 19, 2016 at 4:25 PM, Jonathan Wilkes jancsika@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.
I agree that Pd's code is very sparsely commented. I imagine Miller keeps another copy with all the comments. I was mystified by UNITBIT32 as well. I tried to explain it in the comments of my sqosc~ object at https://sourceforge.net/p/pure-data/svn/HEAD/tree/trunk/externals/mrpeach/sq...
Martin
On Wed, Oct 19, 2016 at 4:25 PM, Jonathan Wilkes jancsika@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.
And if we're not talking hot code then the answer of which answer arrived in less time should (at least in large part) dictate which method to use, shouldn't it?
-Jonathan
Katja
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@gmail.com An: "Jonathan Wilkes" jancsika@yahoo.com Cc: "pd-list@lists.iem.at" pd-list@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@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@gmail.com To: "pd-list@lists.iem.at" pd-list@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@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@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
Pd-list@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
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@gmx.at An: katja katjavetter@gmail.com Cc: "pd-list@lists.iem.at" pd-list@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@gmail.com An: "Jonathan Wilkes" jancsika@yahoo.com Cc: "pd-list@lists.iem.at" pd-list@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@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@gmail.com To: "pd-list@lists.iem.at" pd-list@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@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@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
Pd-list@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
Pd-list@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
On 10/19/2016 11:40 PM, Christof Ressi wrote:
I guess, you can't change the cutoff frequency so often that the CPU cost will actually matter. But maybe I'm wrong...
well, the perform routing will be called *at most* for each block (which can be 1 sample when running with [block~ 1])
the cutoff frequency will be changed whenever you send a message, which can be much ... more often, [until] somebody fixes your patch.
fmar IOhannes
Good point. :-) so let's put it this way: the difference in cpu cost is neglectible for well designed patches. Again, i might be wrong.
Gesendet: Mittwoch, 19. Oktober 2016 um 23:54 Uhr Von: "IOhannes m zmölnig" zmoelnig@iem.at An: pd-list@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 10/19/2016 11:40 PM, Christof Ressi wrote:
I guess, you can't change the cutoff frequency so often that the CPU cost will actually matter. But maybe I'm wrong...
well, the perform routing will be called *at most* for each block (which can be 1 sample when running with [block~ 1])
the cutoff frequency will be changed whenever you send a message, which can be much ... more often, [until] somebody fixes your patch.
fmar IOhannes
Pd-list@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
Hi,
On 19/10/16 15:25, Jonathan Wilkes via Pd-list wrote:
When implemented in C, which approach takes the least amount of time to read, reason about, and fully comprehend?
Agreed, as changing filter frequency at message rate is probably a relatively cold code path, I vote for something like:
double w = 2 * pi * fmin(fmax(fc / SR, 0), 0.5); double k = (1 - sin(w)) / cos(w);
Maybe t_sample would be better than double, but then you'd need C99 type-generic maths or other similar macro tricks (or C++) to call sin or sinf as appropriate.
From: katja katjavetter@gmail.com 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).
Yes that's sensible, due to the odd symmetry.
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.
I went up to 9th degree polynomial in my graphical analysis below.
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.
Yes a fun challenge.
Here's a graph, it shows that polynomial approximations do quite poorly for this problem, even when increasing the degree: https://mathr.co.uk/misc/2016-10-21_filter_coefficient_calculation_accuracy.... (logarithmic frequency scale, if SR is 48kHz it's from ~11Hz to 24kHz)
I used this Haskell to generate the data file for curve fitting https://mathr.co.uk/misc/2016-10-21_filter_coefficient_calculation_accuracy.... and this GNUPlot to fit the curves and plot the graph: https://mathr.co.uk/misc/2016-10-21_filter_coefficient_calculation_accuracy....
Having said that the polynomial approximations do poorly, they're probably accurate enough in practice, I suspect more serious problems would be likely to occur from using single precision in the recursive feedback path in the filter, see for example: https://lists.puredata.info/pipermail/pd-list/2010-08/082104.html
Today pd is even deployed on embedded devices so the frugal coding approach is still relevant. After 20 years.
In my tests last year and revisited in July I found a 9th degreee odd polynomial was both more accurate and more efficient than the tabfudge stuff for linear interpolated cosine table lookup on every single architecture I have available to me (from pentium-4 to raspberry pi 3 via amd64), when compiled with -march=native on the target machine: https://mathr.co.uk/blog/2016-07-20_approximating_cosine_update.html https://mathr.co.uk/blog/2015-04-21_approximating_cosine.html
Claude
-
Christof Ressi -
Claude Heiland-Allen -
cyrille henry -
IOhannes m zmölnig -
Jonathan Wilkes -
katja -
Martin Peach