🚀 go-pugleaf

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

In article <1989Sep28.161516.10353@rpi.edu> wrf@mab.ecse.rpi.edu (Wm Randolph Franklin) writes:

   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.


wrf missed the point entirel, i think.  what charlie was wanting is a
way to compute the amount of distortion present in a digitized sine
wave. 

happily, what he wants is relatively easy to provide.  the fourier
transform does exactly the right thing.  in the discrete time case,
you will have to either sample for a long time, or know the freqency
of the original sine wave very accurately.

so... the answer to the question is to computer the fourier transform,
eliminate the coefficient corresponding to the sine wave you are
interested in and compute the total power in the rest of the spectrum.

--
ted@nmsu.edu
			remember, when extensions and subsets are outlawed,
			only outlaws will have extensions or subsets

Wm Randolph Franklin writes:

> 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.

He goes on to discuss fast ways to compute sines of angles.  First of all, I
don't know how Charlie is computing THD/Noise in a digitized sine wave, but
it doesn't sound like he is trying to generate a sine wave.  Maybe I'm all
wet.  Maybe he knows the "exact" frequency and amplitude and all he needs to
do is quickly calculate the sine wave.  If it has a constant freq/amp, then
it is far easier to simply step through a table.  In fact, even if the
frequency is different (than the one in your lookup table), you just change
your step size.  Finally, even if the amplitude is different (than the one in
your lookup table), you can easily scale the values read from the table by:

  1 - A quick multiplication (if that operation is quick enough), or

  2 - Use a sine function rule which states:

           sin(a-b) + sin(a+b) = f(b) * sin(a)

      where a is the current angle you are trying to compute the sine
      value of and b is a constant "delta" (hence f(b) is also a constant).
      Sorry, I don't know what the formula for f(b) is off the top of my
      head.  It can be derived quite easily, however.  If there is any
      interest, I can look it up...

=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=

  Edward McClanahan
  Hewlett Packard Company
  Mail Stop 47UE              -or-     edwardm%hpda@hplabs.hp.com
  19447 Pruneridge Avenue
  Cupertino, CA  95014                 Phone: (408)447-5651

> charlie@oakhill.UUCP    Charlie Thompson at Motorola Inc. Austin, Tx

> 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.  The fitting algorithm with adapt to phase frequency, and
> amplitude and create a minimum RMS error replica of the digitzed
> sinewave.  The RMS error thus indicates error in the original digitzer.

A good paper on what you are looking for is in an Product Note 5180A-2
"Dynamic Performance Testing of A to D Converters".  The code is in
appendix IV.  This app note includes some source code for what you are
looking for, but it's not in C.  Embarrasingly enough it's in HPL.

This is a good A/D testing technique.  The "effective bits" number is a
good way to compare A/D converters.  Unfortunately there is not good
agreement on what "effective bits" means, but when measured with the
technique you describe it is quite enlightening.

Different A/D testing techniques reveal different types of problems with
A/D's and S/H's.  Anyone who needs to test A/D's will find a treasure
trove of techniques in 5180A-2.

If you find or write the C code, please post or mail it to me.

Tom Anderson       Hewlett-Packard     Signal Analysis Division
toma@hpsad.hp.com              "It's only hardware"
Opinions expressed are my own and not HP's.

> 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.  The fitting algorithm with adapt to phase frequency, and
> amplitude and create a minimum RMS error replica of the digitzed
> sinewave.  The RMS error thus indicates error in the original digitzer.

The requirement can be restated as: Find A, B, C, and F to minimize:

    sum [ (signal(i)-fit(i))**2 ]

where fit(i) = A + B*sin(2*pi*F*time(i)) + C*cos(2*pi*F*time(i))

If F is fixed, a standard linear algebra package should let you find A, B and
C easily.  Often in A/D testing, you know exactly the ratio of sample rate to
signal frequency, and fiddling with F gives an optimistic result.  If you
really want to fit F as well, try getting the approximate frequency with an
FFT, then search for a minimum of distortion versus F.

This method has several advantages over an FFT:

1.  All of the data points are weighted equally.  Most FFTs need to be
    windowed.

2.  The data need not be spaced evenly, as long as you know the time of
    each sample.

3.  You can get reasonable answers with only a fraction of a cycle of
    the input frequency.

I have an alpha version of a stand-alone program that takes ASCII values (or
time-value pairs) and gives distortion numbers, which is available on
request.

Bob Jewett   jewett@hplabs

In article <1989Sep28.161516.10353@rpi.edu> wrf@mab.ecse.rpi.edu (Wm Randolph Franklin) writes:

[many lines deleted]

>You can   work out  how many  Taylor  terms  it would  take  to get that
>accuracy.  Editorial: Taylor series are for wimps.

I'm under the impression that in any given range which doesn't have
singularity points, Chebyshev polynomials (of the first kind) can be
used to "telescope" Taylor series of degree N down to degree N-1
while retaining the same maximum error (but this error may pop up in
places where it wasn't present in the original--kinda like squeezing
a balloon).  The difference between N and N-1 is only ONE term
(although for some common trig functions which have only even- or odd-
powered terms, it may LOOK like two?); which seems rather "wimpy" to me.

I can see that if I'm going to calculate a particular series a jillion
times, one less calculation would be nice, but I don't see that as
being particularly significant otherwise.  Do other Chebyshev
polynomials provide even better telescoping?

[more lines deleted]

--Myron.
--
Myron A. Calhoun, PhD EE, W0PBV, (913) 532-6350 (work), 539-4448 (home).
INTERNET: mac@ksuvax1.cis.ksu.edu
BITNET:   mac@ksuvax1.bitnet
UUCP:  ...{rutgers, texbell}!ksuvax1!harry!mac

[mac@ksuvax1.cis.ksu.edu posts a misunderstanding about telescoping
power series with Chebychev polynomials.]

Actually, if one could always reduce the degree of a polynomial
approximation by one, retaining the same max error, then repeating the
process would give us a zeroth degree approximation with the same max
error as the original Nth degree polynomial provided.
Doesn't seem reasonable.

The method looks more like this:
Generate a Taylor series giving significantly better accuracy ofapproximation
than is needed over the interval of interest, say \epsilon / 2.
Express this (N-th degree) polynomial as a weighted sum
of Chebychev polynomials:
\sum_{k=0}^N \alpha_k T_k ( x).
This step is just a change of basis in the vector space of N-th degree
polynomials.
Now choose M such that \sum_{k=M+1}^N \| \alpha_k \| <= \epsilon / 2,
i.e. the tail of the series is small enough.
Then we use the M-th degree polynomial \sum_{k=0}^M \alpha_k T_k ( x)
as our approximating polynomial, and the error will be less
than or equal to \epsilon / 2 + \epsilon / 2 = \epsilon.

Generally M will be much smaller than N; one is trading off
accurate approximation of the function derivatives near one point
for accurate approximation of function values over the whole interval.


Vincent Broman,  code 632, Naval Ocean Systems Center, San Diego, CA 92152, USA
Phone: +1 619 553 1641    Internet: broman@nosc.mil   Uucp: sdcsvax!nosc!broman

You're right that my calling Taylor series wimps  is excessive.  This is
especially  so for  sin(x) which  is  one of  the functions  that Taylor
series are best  suited for  since sin  is smooth  and defined  over the
whole complex plane with  no singularities, and  the Taylor coefficients
shrink fast.

However,  even  here Chebyshev polynomials are  better, and the relative
difference increases as the order increases.

I scaled sin to that [0,Pi] was mapped to [-1,1]  and expanded about the
new  0 (for Taylor)   or over the interval  (for   Chebyshev).  This  is
probably the interval you would want to approximate.

N	abs(Coef of x^n) in Taylor	 	
                   abs(Coef of T(N,x)) in Cheb
0	.7        .6   
1	.55	  .5   
2	.22	  .1   
3	.06	  .01  
4	.01	  .001 
5	.002	  .0001
6	.0002	  7e-6 
7	.00002	  3e-7 
8	3e-6	  2e-8 
9	2e-7	  9e-10
10	2e-8	  3e-11
11	1e-9	  1e-12
12	8e-11	  4e-14
13	5e-12     1e-15
14	3e-13	  3e-17
15	1e-14	  9e-19

Note: 

1. Truncating the Chebyshev series when it is still expressed  in T(n,x)
is equivalent (almost) to truncating the Taylor series expressed in x^n.
In either case, the error from stopping at any term is usually less than
the first unused coefficient.  (Picky note to Chebyshev hackers: this is
why  it is more  useful to  look  at the   coefficients of T rather than
convert to a power basis.)

2. Although I only retyped 1 or 2  digits for this  note, the calcs were
done to 20 significant digits to avoid roundoff.

3. The Taylor to the 15th  degree is  between the 12th  and  13th degree
Chebyshev -- which is more than one term.


Well,  maybe  you think that   the interval   favored Chebyshev. So try
another interval [-Pi,Pi], scaled to [-1,1].  Note the function is odd.

N Taylor   Cheb
1  3      1     
3  5      1.6   
5  3      .2    
7  .6     .01   
9  .08    .0005 
11 .007   .00001
13 .0004  2e-7  
15 .00002 2e-9  

The Taylor to x^15 worse than the Chebyshev to x^11.

----------------

In all the  above I'm using  about the best   possible function for  the
Taylor.   In fact Chebyshev   can do a  wide class  of  functions Taylor
can't:

- functions with a complex singularity nearby, such as 1/(x^2+.1).

-  functions  which are  only differentiable  a  small number  of times.
Taylor requires N for N-th degree, Chebyshev requires only 1st order.

Admittedly Chebyshev will converge much more slowly in these cases.  But
it WILL converge, and Taylor WON'T.

----------------

(Now  I'm really getting on my  hobby-horse).   There's  something often
much better than  Chebyshev --  a Pade-Chebyshev expansion,  which  is a
Chebyshev  converted to  a rational  fraction.  Its  relative betterness
also increases with N.  However it requires one divide to evaluate.  The
problem is that Pade expressions are not  in the textbooks yet, and even
if they were are more complicated to calculate.  Also  I can't get Maple
to do it right now.

-- 
						   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

In article <TED.89Sep28182735@aigyptos.nmsu.edu> ted@nmsu.edu (Ted Dunning) writes:
>   >used to compute THD/Noise in a digitized sine wave.
>
>way to compute the amount of distortion present in a digitized sine
>
>so... the answer to the question is to computer the fourier transform,
>eliminate the coefficient corresponding to the sine wave you are
>interested in and compute the total power in the rest of the spectrum.
>
Yes, this appears to be a feasible solution to the problem, except for
artifacts introduced by the DFT/FFT algorithm.  The DFT (read FFT) is
able to work because it assumes that the time series that is input to 
the DFT is periodic.  In other words, it assumes that after all the 
input samples have been scooted into the DFT buckets for processing, that
an identical waveform will occur in the same N time samples immediately 
following.  If your sine wave has a period that exactly coincides
with the number of samples in your DFT (or an integer division thereof)
then the output spectrum will have a single positive (and negative) 
frequency component at the frequency of the input sine wave.

If the sine wave does not have an integer number of periods in the 
input time sequence to the DFT, then a smearing will appear in the 
spectral characteristic.  This smearing, or leakage, is due to a discontinuity
in the edges of the input time sequence.  Whereas a sine wave with integer
periods will be continuous when the input N point sequence is lined up
end to end, having a sine wave with a period that does not allow
for an integer number of repetitions of the sine wave in the input
buckets of the DFT will have an abrupt discontinuity at the abbutment point.

The DFT calculates a spectrum that attempts to recreate this discontinuity.

Therefore, the use of windows may be called for.  Windowing the input sequence
will squush the discontinuities at the ends of the input time sequence.  The 
standard reference for windowing is:

	Harris F.J. 1978, On the use of Windows for Harmonic Analysis 
 	with the DFT, Proc IEEE, 66 (January)

The idea with the Chebyschev sine wave computation may also open a possibility
toward determining the THD.  If a good clean sine wave, of exactly the
same period and phase of the signal in question may be generated, then the
difference may be taken, and the spectrum of the difference transformed
to obtain the signal of which to determine the power of the harmonic    
distortion.  

Good Luck

Greg Harrison
My opinions are not intended to reflect those of GE.

In <3811@deimos.cis.ksu.edu> mac@harris.cis.ksu.edu (Myron A. Calhoun) writes:
>I'm under the impression that in any given range which doesn't have
>singularity points, Chebyshev polynomials (of the first kind) can be
>used to "telescope" Taylor series of degree N down to degree N-1
>while retaining the same maximum error (but this error may pop up in
>places where it wasn't present in the original--kinda like squeezing
>a balloon).

Well, it was extreme of me to call Taylor series wimpy, especially since
sin  is  about   the best example   of a  Taylor  --   it's  entire  (no
singularities)  and the   coefficients decrease fast.    Also,  a Taylor
series is often by far the easiest approximation to a function.

However, even with sin, Chebyshev have an advantage of more than 1 term,
and the advantage increases with the degree.   Consider sin(x) from 0 to
pi/2 -- the obvious interval for a  sin implementation.   I expanded the
Taylor series about pi/4, and  did the Chebyshev  over [0,pi/2], both to
20 significant  digits.  These coefs are  the approx  max  error at each
stage since I also renormalized x to be in [-1,1].

   n abs(Taylor coef) abs(Cheb coef (i.e. coef of T(n,x)))
   5 2e-3        7e-6
   10 2e-8	 3e-11
   15 1e-14	 9e-19

The 15th degree Taylor has an accuracy  between  the 12 and  13th degree
Chebyshev -  more  than 1   degree  improvement.  For  other series   or
intervals, it can be more.  Consider sin(x) from -pi  to pi.  I expanded
sin(x)/x  about 0.  Note that  it's an  even function so there's  no 5th
degree.

   n Taylor Chebyshev
   4 3      .2
   10 7e-3 1e-5
   14 2e-5 2e-9

The 14th degree Taylor is between an 8th and 10th degree Chebyshev.

Next,  your restriction to non singular  functions is quite strict since
many smooth, infinitely differentiable, functions, such as 1/(x^2+0.01),
have a radius  of convergence for  the  Taylor  of  0.1  because of  the
complex pole  at 0+.1i.  Also,  an N-th degree  Taylor requires that the
function be differentiable  N   times, whereas  Chebyshev  requires only
once.  

In both  of these cases,  the  Chebyshev series converges slowly, but it
converges over the whole   interval,  whereas  Taylor converges  only in
(-0.1, 0.1) in the first, and not at all in the second.

Actually, if you allow one division in your approximation, you can use a
rational or Pade form of either the Taylor or Chebyshev, which  for less
well behaved functions vastly improves either  approx. They're much more
complex to calculate, but can improve on the Taylor even for exp(x).

Summary: Taylor is certainly ok for sin(x), but Chebyshev is faster even
here, and does a much wider class of functions.  However,  it is messier
to calculate.


-- 
						   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