Sine Lookup Table with Linear Interpolation

The simple method would be to compute a zillion-point FFT. Since the function is really truly periodic, the FFT gives you the right answer for the continuous-time Fourier transform, provided you have enough points that you can see all the relevant harmonics.

Cheers

Phil Hobbs

Oops, should be: cos(x) ~= 1 - x^2/2 sin(x) ~= x (- x^3/6) Mixed up your odd/even powers :)

Of course, the cubic is less important than the quadratic, so sin(x) ~= x is sufficient for most purposes.

I suppose you could potentially synthesize the functions using low power approximations (like these) and interpolating between them over modulo pi/4 segments. Probably wouldn't save many clock cycles relative to a CORDIC, higher order polynomial, or other, for a given accuracy, so the standard space-time tradeoff (LUT vs. computation) isn't changed.

Tim

harmonic

just

component

null

Why can't some harmonic component be added to the table points to improve downstream distortion? Of course, you'd need a full 360 degree table, not a folded one.

filter,

table

For everything, including the DAC, but especially downstream stuff. The amps and stuff may have 100x more distortion than the lookup table does.

All you need now is a digitizer or spectrum analyzer that has PPM distortion!

We've bumped into this problem, trying to characterize some of our ARBs. The best thing is to use a passive lowpass or notch filter (with very good Ls and Cs!) to remove the fundamental and analyze what's left.

hehe. yes. But it gets the point across, either way. Detail I didn't care a lot about at the moment.

Yup.

Jon

In comp.dsp rickman wrote: (snip)

The cos(M)sin(L) is one table, with some bits from each.

Well, when you do the interpolation, you first need to know the spacing between points in the sine table. That spacing is proportional to cos(x). So, the interpolation table is indexed by the high four bits of M and the four bits of L.

(snip on getting to 90 degrees, or 2^n in binary.

(snip)

I am not sure yet how they do that. I did notice that when adding the interpolation value they propagate the carry, so the output can go to the full 1.0000 instead of 0.99999. But it will do that before 90 degrees. (snip)

So, the one referenced uses three 256x4 ROMs to get 12 bits of sin(M) from 8 bits of M, and then a fourth 128x8 ROM to get five bits to add to the low bits of the 12 bit sin(M).

-- glen

i think about hardware as a technology, not so much an app. i think that apps can have either hardware or software realizations.

in *software* when using a LUT *and* linear interpolation for sinusoid generation, you are interpolating between x[n] and x[n+1] where n is from the upper bits of the word that comes outa the phase accumulator:

#define PHASE_BITS 10 #define PHASE_MASK 0x001FFFFF // 2^21 - 1 #define FRAC_BITS 11 #define FRAC_MASK 0x000007FF // 2^11 - 1 #define ROUNDING_OFFSET 0x000400

long phase, phase_increment, int_part, frac_part;

long y, x[1025];

...

phase = phase + phase_increment; phase = phase & PHASE_MASK;

int_part = phase >> FRAC_BITS; frac_part = phase & FRAC_MASK;

y = x[int_part]; y = y + (((x[int_part+1]-y)*frac_part + ROUNDING_OFFSET)>>FRAC_BITS);

now if the lookup table is 1025 long with the last point a copy of the first (that is x[1024] = x[0]), then you need not mask int_part+1. x[int_part+1] always exists.

uh, because 2^10 = 1024?

that came from your numbers: 8-bit address for a single quadrant and there are 4 quadrants, 2 bits for the quadrant.

i'm agnostic about the specs for your numbers. you are saying that 8 bits per quadrant is enough and i am stipulating to that.

why wouldn't you use those extra bits in the linear interpolation if they are part of the phase word? not doing so makes no sense to me at all. if you have more bits of precision in the fractional part of the phase

not in software. i realize that in hardware and you're worried about real estate, you might want only one quadrant in ROM, but in software you would have an entire cycle of the waveform (plus one repeated point at the end for the linear interpolation).

and if you *don't* have the size constraint in your hardware target and you can implement a 1024 point LUT as easily as a 256 point LUT, then why not? so you don't have to fiddle with reflecting the single quadrant around a couple of different ways.

i read this over a couple times and cannot grok it quite. need to be more explicit (mathematically, or with C or pseudocode) with what your operations are.

not with defining the points of the table. you do that in MATLAB or in C or something. your hardware gets its LUT spec from something you create with a math program.

i'm just thinking that if you were using LUT to generate a sinusoid that is a *signal* (not some parameter that you use to calculate other stuff), then i think you want to minimize the mean square error (which is the power of the noise relative to the power of the sinusoid). so the LUT values might not be exactly the same as the sine function evaluated at those points.

what are "lines"? not quite following you.

if you want to define your LUT to minimize the mean square error, assuming linear interpolation between LUT entries, that can be a little complicated, but i don't think harder than what i remember doing in grad school. want me to show you? you get a big system of linear equations to get the points. if you want to minimize the maximum error (again assuming linear interpolation), then it's that fitting a straight line to that little snippet of quadratic curve bit that i mentioned.

sin(t) = t + ... small terms when t is small

cos(t) = 1 - (t^2)/2 + ... small terms when t is small

but they beat you by a few hundred years.

:-)

If you are implementing in an FPGA with no multipliers, then CORDIC is the first thing you should be looking at. Do a web search on "CORDIC" for references. I know that's not at all obvious to you, but it works to try. Or you could keep coming back here for more spoonfeeding.

Eric Jacobsen Anchor Hill Communications

Problem with that..could see only front and back cover. Different source perhaps?

(snip, I wrote)

I just tried it, copy and paste from the above link, and it worked.

You might be sure to use Adobe reader if that matters.

Also, it is page 273 in the document, page 282 in the PDF.

(snip)

Once upon a time, about 30 years ago, i fiddled with routines to multiply a millions-of-digits number by another millions-of-digits number. Naturally, i used the FFT and convolution (and bit-reordering) to do this. As the number of digits grew, i needed to use a larger N (in 2^N samples), so different algorithms were needed in different ranges to maximize speed. So i had to try a goodly number of published, standard, accepted, used routines. ALL algorithms for N larger than 5 worked but i was able to speed them up anywhere from 10% to 30% by using "tricks" which were the equivalent of "folding" around 22.5 degrees. Start with a full 360 degree method. "Fold" in half: sin(-x)=-sin(x) or 180 degrees. Half again on the imaginary axis for what you are doing for 90 degrees: sin(90-x) = sin(90)*cos(x)-cos(90)*sin(x) = -cos(x) exact. Half again: sin(45)=cos(45)=2^(3/2)~~.707 carry out as far as you want; sin/cos to 45 degrees. sin(x=0..45) direct, sin(x=45..90) =cos(45-x) accuracy depending on value of sin/cos(45). That is prolly as far as you want to go with a FPGA, because (as you indicated) you do not want to attempt multiplies in the FPGA. That last step of "folding" might give better results or better FPGA / LUT calcs.

Calculating actual values for the LUT for maximum accuracy: the trig series x-(x^3/3!)+(x^5/5!) etc should give better results from 0..90 than the approx on page 51. Do not hesitate to try alternate schemes: sin(2x) = 2*sin(x)*cos(x) and sin(0.5*x)=+/- sqrt((1-cos(x))/2) where that sign depends on the quadrant of x.

?? CORDIC requires even more multiplies

?? CORDIC requires even more multiplies

But you can muliply two short lookup tables to simulate a longer lookup table using he identity:

sin(a+b) = sin(a)*cos(b) + cos(a)*sin(b)

so you do a coarse lookup table with sin(a) (and you can get cos(a) from it ti if a is some integer factor of 90 degrees)

and do two fine lookup tables for sin(b) and cos(b) eg: you can do 16 bits of angle at 16 bit resoluton in 512 bytes of table and ~80 lines of assembler code (but the computation needs two multiplies) the first two bits determine the quadrant. the next 7 bits are on a coarse lookup table 128 entries x 16 bits is 256 bytes next bit defines the sign of b (if set b is negative,add one to a) and the final 6 are the rest of b from a 64 entry table you can get sin(|b|) and another 64 for cos(|b|) treat b=0 as a special case - have entries for 1..64 in the table correct sin(|b|) for negative b a few branches and bit-shifts, two multiplies, one add and you have your value.

CORDIC requires shift and add. That's why it was used in early scientific calculators- very little hardware was required.

I have an HP 339A you can borrow. ;)

Cheers

Phil Hobbs

That's for audio. Who cares about distortion for audio? And it gets nowhere near PPM turf.

rickman schrieb:

Hello,

you may use the well known formula: sin(a+b) = sin(a)cos(b) + cos(a)sin(b)

It looks like you will need four tables, one coarse sin(), one coarse cos(), one fine sin() and one fine cos() table, but you only need the coarse sin() for one quadrant and the fine sin() and cos() tables. The coarse sin() table for one quadrant will also deliver the values for cos()

Bye

Just one note, since you're doing this in an FPGA. If your look up table is a dual-port block RAM, then you can look up the the slope simultaneously rather than calculate it. sin(x)/dx = cos(x) = sin(x +

It works for me. Downloaded & opened with Adobe Acrobat X Ver.10.1.6