charlie@oakhill.UUCP (Charlie Thompson) (09/28/89)
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. Thanks, Charlie Thompson
wrf@mab.ecse.rpi.edu (Wm Randolph Franklin) (09/28/89)
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
ted@nmsu.edu (Ted Dunning) (09/29/89)
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
edwardm@hpcuhc.HP.COM (Edward McClanahan) (09/29/89)
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
toma@hpsad.HP.COM (Tom Anderson) (09/30/89)
> 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.
jewett@hpl-opus.HP.COM (Bob Jewett) (09/30/89)
> 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
harrison@sunny.DAB.GE.COM (Gregory Harrison) (10/06/89)
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.