Hello everyone
I am trying to modify fir1.m to create FIRs using the correct sum of sines or cosines. I ended up making a new .m file, but when I got to the "%! assert()" lines, I stumbled across some numerical issues. The .m file was built after a wxMaxima file, which was used to generate the default results, against which assert() could check, and I noticed that both sin(), and cos(), seem to not produce (anti-)symmetrical results. The differences can be more than 1e-15. "format long" is enabled, but it's the same. Example, lowpass, type I, 33rd order, wc=0.37, in wxMaxima and Octave (in this order): 0.006248855546642409" " -0.0151891541716538" " -0.02000749470712845" " 3.703551396872747*10^-4" " 0.02352639910729612" " 0.01987710435709456" " -0.0107156702183056" " -0.03346910762203006" " -0.01647493965701502" " 0.02756343124616879" " 0.0468059024239567" " 0.006350821104335874" " -0.06144311077699433" " -0.07272786646455694" " 0.02972318687964276" " 0.2090466916579446" " 0.3495187814185787" " 0.3495187814185787" " 0.2090466916579446" " 0.02972318687964276" " -0.07272786646455694" " -0.06144311077699433" " 0.006350821104335874" " 0.0468059024239567" " 0.02756343124616879" " -0.01647493965701502" " -0.03346910762203006" " -0.0107156702183056" " 0.01987710435709456" " 0.02352639910729612" " 3.703551396872747*10^-4" " -0.02000749470712845" " -0.0151891541716538" " 0.006248855546642409 6.24885554664247e-03 -1.51891541716537e-02 -2.00074947071285e-02 3.70355139687233e-04 2.35263991072961e-02 1.98771043570946e-02 -1.07156702183055e-02 -3.34691076220301e-02 -1.64749396570151e-02 2.75634312461687e-02 4.68059024239567e-02 6.35082110433597e-03 -6.14431107769943e-02 -7.27278664645570e-02 2.97231868796425e-02 2.09046691657944e-01 3.49518781418579e-01 3.49518781418577e-01 2.09046691657945e-01 2.97231868796429e-02 -7.27278664645570e-02 -6.14431107769943e-02 6.35082110433577e-03 4.68059024239567e-02 2.75634312461688e-02 -1.64749396570150e-02 -3.34691076220301e-02 -1.07156702183056e-02 1.98771043570945e-02 2.35263991072961e-02 3.70355139687233e-04 -2.00074947071284e-02 -1.51891541716538e-02 6.24885554664247e-03 The simplest generator for this would be (sinc() works the same): N=33; M=N/2; n=-M:1:M; h=sin(0.37*pi*n)./(pi*n); h(ceil((N+1)/2)=0.37; Note that the code in Octave directly follows the one in wxMaxima, and both use the range -M:1:M. There is no cheating in wxMaxima, so sin(-1)=-sin(1), while in Octave it doesn't seem to be so. What's more, I am not using big floats in wxMaxima, only default numbers, float() when necessary. The assert() that I was aiming for was 1e-15. In the case above, it fails at h(18): Location | Observed | Expected | Reason (18) 0.34952 0.34952 Abs err 1.6653e-15 exceeds tol 1e-15 Is this a known behaviour? Should I just relax the tolerance to 1e-14, or less? Or am I doing something wrong? I'll attach the .m script as it is, so far. Vlad fir.m (11K) Download Attachment |
I have just seen a previous question in the mailing lists, "If pi is
so accurate why it's not producing that accurate result". There probably are more of the same. Just to clarufy (if there's a misunderstanding), I am not questioning the inherent floating point errors, but their consistency. In Octave, sin(-x) doesn't seem to be equal to -sin(x). This is my main point. As you can see, the the order of the FIR is odd, so it has two midpoints of equal values coming from sin(x)/x and sin(-x)/(-x), thus, for all intents and purposes, sin(x)/x, yet they are not equal. That they differ in Octave and wxMaxima, so be it, to each their own (numerical evaluation), but why the two are different whithin Octave is my concern. Vlad |
In reply to this post by vlad ionescu
On Mon, Mar 26, 2018 at 16:01:04 +0000, vlad ionescu wrote:
> The simplest generator for this would be (sinc() works the same): > > N=33; M=N/2; n=-M:1:M; h=sin(0.37*pi*n)./(pi*n); h(ceil((N+1)/2)=0.37; > > Note that the code in Octave directly follows the one in wxMaxima, and > both use the range -M:1:M. There is no cheating in wxMaxima, so > sin(-1)=-sin(1), while in Octave it doesn't seem to be so. Setting the fir function aside, I don't see any problems with this particular example >> N = 33; M = N/2; n = (-M:M)'; h = sin (0.37 * pi * n) ./ (pi * n); >> assert (h, flipud (h)) >> assert (h, x, 1e-15) It doesn't look to me like the errors are in the sin function. FWIW the h I obtain above is not equal to the h I get from the attached fir function. The error is worst at index 18, the same as the error you are seeing. -- mike signature.asc (849 bytes) Download Attachment |
On Tue, Mar 27, 2018 at 11:29:54 -0700, Mike Miller wrote:
> It doesn't look to me like the errors are in the sin function. In fact, I think what you want to look at are the differences between the following v1 and v2 N = 33; M = N/2; n = -M:M; v1 = sin (0.37 * pi * n) ./ (pi * n); v2 = sin (0.37 * pi * [n]) ./ (pi * [n]); I think what is hurting your precision here is the lazy evaluation of ranges in Octave. -- mike signature.asc (849 bytes) Download Attachment |
I ran across a fairly recent (2008) paper on some of the "dark arts" of
finding good approximations to mathematical functions. It does in part look at sin(x). The paper has a preferred citation: Florent De Dinechin, Christoph Lauter. Optimizing polynomials for floating-point implementation. 12 pages. 2008. <ensl-00260563> I believe the URL is: https://hal-ens-lyon.archives-ouvertes.fr/ensl-00260563 Gord |
________________________________ From: Gordon Haverland <[hidden email]> To: [hidden email] Sent: Wednesday, March 28, 2018 1:53 AM Subject: Re: Numerical issues with sin()/cos() I ran across a fairly recent (2008) paper on some of the "dark arts" of finding good approximations to mathematical functions. It does in part look at sin(x). The paper has a preferred citation: Florent De Dinechin, Christoph Lauter. Optimizing polynomials for floating-point implementation. 12 pages. 2008. <ensl-00260563> I believe the URL is: https://hal-ens-lyon.archives-ouvertes.fr/ensl-00260563 Gord ----------------------------------------------------- Since Intel 80387 (see https://en.wikipedia.org/wiki/X86_instruction_listings -> https://en.wikipedia.org/wiki/X86_instruction_listings#Added_with_80387 ) sine and cosine functions are part of the instruction set, so if one wants fast implementation, he's limited by what his x86 (if it's the case) provides. --Sergei. |
In reply to this post by vlad ionescu
Hello
I am sorry, despite being able to see all the (current) 4 replies in the gnu.org list, Gmail still didn't bother to let me know, so I am replying to my own sent mail. > v1 = sin (0.37 * pi * n) ./ (pi * n); > v2 = sin (0.37 * pi * [n]) ./ (pi * [n]); > > I think what is hurting your precision here is the lazy evaluation of > ranges in Octave. This is the first time I see this, so I did some search. Correct me if I am wrong, but, apparently, it's about the way Octave stores the range. If it's without square braces, it uses less storage, which probably means it uses floats instead of doubles, which might explain the precision (in fact, I am pleasantly surprised by the accuracy, now), but that still leaves the question of why sin(-x) != -sin(x), though I guess it's the floating point representation of the numbers, as they go up, or down, different fractions (just guessing), but then... > FWIW the h I obtain above is not equal to the h I get from the attached > fir function. The error is worst at index 18, the same as the error you > are seeing. Could it be that different CPU architectures calculate differently the floats? I can't imagine why you would get different results. Did you use fir(33,[0 0.37])? Still, the eror is, again, at index 18, so maybe it really is about the floating point representation. At any rate, as soon as I made my range as n=[...], everything went spot on. Thank you. > I ran across a fairly recent (2008) paper on some of the "dark arts" of > finding good approximations to mathematical functions. It does in part > look at sin(x). Well, for my part, I'm fine with the way Octave does it, thank you for the suggestion. Vlad |
On Wed, Mar 28, 2018 at 06:46:47 +0000, vlad ionescu wrote:
> > I think what is hurting your precision here is the lazy evaluation of > > ranges in Octave. > > This is the first time I see this, so I did some search. Correct me if > I am wrong, but, apparently, it's about the way Octave stores the > range. If it's without square braces, it uses less storage, which > probably means it uses floats instead of doubles Actually it's that Octave stores ranges as a (start, step, stop) triplet. Let's say the value n is the range -16.5:1:16.5 as in your example. Octave stores that as (start=-16.5, step=1, stop=16.5). The expression 0.37*pi*n evaluates to another range, with the original triplet values each multiplied by 0.37*pi. Each of the 34 values of the range has not been evaluated yet. It's only when the range is passed in to the sin function that it is evaluated into an array of discrete doubles. So I think you can see how this lazy evaluation can lead to differences when adding subsequent increments of 0.37*pi to the start value of -0.37*pi*16.5, compared to what you actually wanted. > > FWIW the h I obtain above is not equal to the h I get from the attached > > fir function. The error is worst at index 18, the same as the error you > > are seeing. > > Could it be that different CPU architectures calculate differently the > floats? I can't imagine why you would get different results. Did you > use fir(33,[0 0.37])? Still, the eror is, again, at index 18, so maybe > it really is about the floating point representation. Turns out this was only because I had tested with n=(-16.5:16.5)', where the transpose operator also expands the range into a column vector. -- mike signature.asc (849 bytes) Download Attachment |
Free forum by Nabble | Edit this page |