Article View: comp.dsp
Article #289Re: Sine curve fit algorithms
From: Wm Randolph Franklin
Date: Thu, 28 Sep 1989 16:15
Date: Thu, 28 Sep 1989 16:15
70 lines
4355 bytes
4355 bytes
In <2421@radio.oakhill.UUCP> charlie@oakhill.UUCP (Charlie Thompson) writes:
>I am looking for C source code for sine wave curve fitting algorithms. Such an algorithm can be
>used to compute THD/Noise in a digitized sine wave.
For simplicity a Taylor series would be ok. However, here is a fourth
order Chebyshev series that is much faster for the accuracy. It finds
sin(y) for y from 0 to Pi/2 with a max error of 0.000124. This was
obtained with Maple.
2 3
c2 := .00011389456 + .9962889583 y + .0195161068 y - .2031187635 y
4
+ .02856566868 y
If I plot both sin(y) and c2 you can't tell any difference, so here is
the error plotted:
> plot(sin(y)-c2,y=0..Pi/2);
+ 0.000124
| AAA
| A AA AAAAAAA A
| A A A A A
| A A AA A A
+ 0.0000624 A A AA A
| A A A A A
| A A A A A
| A A A A A
| A A A A A
| A A A A A y
+-0-------------*---+-----------------*+-------------------+--*-------------*-+
0 A A .393 .785 1.18 A 1.57
| A A A A A
|A A A A A
|A A A A A
+A-0.0000624 A A A A
|A A A A A
* A A A A
* AA AA A A
* A AA AAA
* AAA
+ -0.000125
>
You can work out how many Taylor terms it would take to get that
accuracy. Editorial: Taylor series are for wimps.
If the angle is outside [0,1.57] use the obvious reduction formulae.
If you want the formula in degrees, or a different accuracy, tell me and
I'll work it out.
Historical note: IBM in Fortran finds sin(x) by reducing the angle to
[0,1.57], then using a half angle formula to reduce x further to [0,
22.5 degrees in radians], then uses a 2 term Chebyshev. This gives full
single (maybe double) precision accuracy.
Note on Chebyshev series: They give the best maximum error over an
interval, for polynonials of a given degree. However, the error
increases very fast outside the interval, so don't try to enlarge the
interval w/o recalculating the series. Also, if you want reduced
accuracy and increased speed, it is not sufficient to just drop a high
order term or 2. You have to go back to the Chebyshev basis.
--
Wm. Randolph Franklin
Internet: wrf@ecse.rpi.edu (or @cs.rpi.edu) Bitnet: Wrfrankl@Rpitsmts
Telephone: (518) 276-6077; Telex: 6716050 RPI TROU; Fax: (518) 276-6261
Paper: ECSE Dept., 6026 JEC, Rensselaer Polytechnic Inst, Troy NY, 12180
Technical Details
Message-ID:
<1989Sep28.161516.10353@rpi.edu>Path:
arc1-ca.pugleaf.net!utzoo!attcan!utgpu!jarvis.csri.toronto.edu!mailrus!uwm.edu!gem.mps.ohio-state.edu!rpi!wrfReferences:
<2421@radio.oakhill.UUCP>Show Full Headers
Complete Article Headers
Path: arc1-ca.pugleaf.net!utzoo!attcan!utgpu!jarvis.csri.toronto.edu!mailrus!uwm.edu!gem.mps.ohio-state.edu!rpi!wrf From: wrf@mab.ecse.rpi.edu (Wm Randolph Franklin) Newsgroups: comp.dsp,rec.audio Subject: Re: Sine curve fit algorithms Keywords: sine fit algorithms Message-ID: <1989Sep28.161516.10353@rpi.edu> Date: 28 Sep 89 16:15:16 GMT References: <2421@radio.oakhill.UUCP> Organization: Rensselaer Polytechnic Institute, Troy NY Lines: 70