Forums

Sine Lookup Table with Linear Interpolation

Started by rickman March 16, 2013
On Mon, 18 Mar 2013 14:10:54 -0400, Jerry Avins <jya@ieee.org> wrote:

> >>> I see, I assumed you had an MCU (or fpga) with multipliers >>> >>> considered CORDIC? >> >> I should learn that. I know it has been used to good effect in FPGAs. >> Right now I am happy with the LUT, but I would like to learn more about >> CORDIC. Any useful referrences? > >http://www.andraka.com/files/crdcsrvy.pdf should be a good start. Do you >remember Ray Andraka? > >Jerry >-- >Engineering is the art of making what you want from things you can get. >&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;&#2013266095;
Thranx a bunch. But then you are a goto kind of guy. ?-)
On Sun, 17 Mar 2013 18:02:27 -0400, rickman <gnuarm@gmail.com> wrote:

> >>> One thing I learned while doing this is that near 0 degrees the sin() >>> function is linear (we all knew that, right?) but near 90 degrees, the >>> sin() function is essentially quadratic. Who would have thunk it? >>> >> Sounds like you are making excellent improvements on the standard ho-hum >> algorithms; the net result will be superior to anything done out there >> (commercially). >> With the proper offsets, one needs only 22.5 degrees of lookup ("bounce" >> off each multiple of 45 degrees). > >I can't say I follow this. How do I get a 22.5 degree sine table to >expand to 90 degrees? > >As to my improvements, I can see where most implementations don't bother >pushing the limits of the hardware, but I can't imagine no one has done >this before. I'm sure there are all sorts of apps, especially in older >hardware where resources were constrained, where they pushed for every >drop of precision they could get. What about space apps? My >understanding is they *always* push their designs to be the best they >can be.
I agree. 45 degree sine and cosine table will do the job as will a 90 degree sine or cosine table. ?-)
On Sun, 17 Mar 2013 18:44:07 -0400, rickman <gnuarm@gmail.com> wrote:

>On 3/17/2013 6:20 PM, John Larkin wrote: >> >> It would be interesting to distort the table entries to minimize THD. The >> problem there is how to measure very low THD to close the loop. > >It could also be a difficult problem to solve without some sort of >exhaustive search. Each point you fudge affects two curve segments and >each curve segment is affected by two points. So there is some degree >of interaction. > >Anyone know if there is a method to solve such a problem easily? I am >sure that my current approach gets each segment end point to within &#2013266097;1 >lsb and I suppose once you measure THD in some way it would not be an >excessive amount of work to tune the 256 end points in a few passes >through. This sounds like a tree search with pruning. > >I'm assuming there would be no practical closed form solution. But >couldn't the THD be calculated in a simulation rather than having to be >measured on the bench?
The calculation can be done (most spice versions will do this, even excel can be used), but that doesn't mean that bench testing will give the same results. The calculation must be done on the final value series. ?-)
On Sun, 17 Mar 2013 17:22:59 -0700, John Larkin
<jjlarkin@highNOTlandTHIStechnologyPART.com> wrote:

> >>I don't think you are clear on this. The table and linear interp are as >>good as they can get with the exception of one or possibly two lsbs to >>minimize the error in the linear interp. These is no way to correct for >>this by adding a second harmonic, etc. Remember, this is not a pure >>table lookup, it is a two step approach. > >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.
How do you get that? All harmonics would be within the 90 degree table. Harmonic phase shifts is a different matter though, that would take careful characterization of production quantities of the final product, $$$$$. ?-)
In comp.dsp josephkk <joseph_barrett@sbcglobal.net> wrote:

(snip, someone wrote)
>>I'm not sure this would be easier. The LUT an interpolation require a >>reasonable amount of logic. This requires raising X to the 5th power >>and five constant multiplies. My FPGA doesn't have multipliers and this >>may be too much for the poor little chip. I suppose I could do some of >>this with a LUT and linear interpolation... lol
> No multipliers(?); ouch. My Simpson's rule method may fit but it will > need a trapezoidal strip out of a full Wallace tree multiplier. Preferably > three (speed issue, that may not be that big for you). How do you do the > linear interpolation without a multiplier?
The LUT is indexed by the high (4) bits of the input and the low bits to interpolate with. So it is, pretty much, 16 different interpolation tables which, in other words, is multiply by table lookup. -- glen
On Mon, 18 Mar 2013 19:17:53 -0400, Phil Hobbs
<pcdhSpamMeSenseless@electrooptical.net> wrote:

>On 3/18/2013 10:26 AM, John Larkin wrote: >> On Mon, 18 Mar 2013 10:12:45 -0400, Phil Hobbs >> <pcdhSpamMeSenseless@electrooptical.net> wrote: >> >>> On 03/17/2013 08:27 PM, John Larkin wrote: >>>> On Sun, 17 Mar 2013 19:31:38 -0400, Phil Hobbs >>>> <pcdhSpamMeSenseless@electrooptical.net> wrote: >>>> >>>>> On 3/17/2013 6:44 PM, rickman wrote: >>>>>> On 3/17/2013 6:20 PM, John Larkin wrote: >>>>>>> >>>>>>> It would be interesting to distort the table entries to minimize THD. The >>>>>>> problem there is how to measure very low THD to close the loop. >>>>>> >>>>>> It could also be a difficult problem to solve without some sort of >>>>>> exhaustive search. Each point you fudge affects two curve segments and >>>>>> each curve segment is affected by two points. So there is some degree >>>>>> of interaction. >>>>>> >>>>>> Anyone know if there is a method to solve such a problem easily? I am >>>>>> sure that my current approach gets each segment end point to within &#2013266097;1 >>>>>> lsb and I suppose once you measure THD in some way it would not be an >>>>>> excessive amount of work to tune the 256 end points in a few passes >>>>>> through. This sounds like a tree search with pruning. >>>>>> >>>>>> I'm assuming there would be no practical closed form solution. But >>>>>> couldn't the THD be calculated in a simulation rather than having to be >>>>>> measured on the bench? >>>>>> >>>>> >>>>> 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. >>>>> >>>> >>>> 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. >>>> >>>> >>> >>> 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. >> >> > >You're no fun anymore. :-) > >It does have a 0.01% FS range, so you'd expect to get a lot better than > -80 dB. > >But for fast stuff, I'd probably do the same as you--filter the >daylights out of it and amplify the residuals. > >Cheers > >Phil Hobbs
Which is how most old school distortion analyzers work. Pop'tronics did a four part series on one that used auto-tuned parallel-t filters with over 80 db fundamental notch rejection in the late 60s or early 70s. Realy a slick design in a lot of ways, especially for a hobby mag. ?-)
On Thu, 21 Mar 2013 08:02:56 -0700, josephkk <joseph_barrett@sbcglobal.net>
wrote:

>On Mon, 18 Mar 2013 19:17:53 -0400, Phil Hobbs ><pcdhSpamMeSenseless@electrooptical.net> wrote: > >>On 3/18/2013 10:26 AM, John Larkin wrote: >>> On Mon, 18 Mar 2013 10:12:45 -0400, Phil Hobbs >>> <pcdhSpamMeSenseless@electrooptical.net> wrote: >>> >>>> On 03/17/2013 08:27 PM, John Larkin wrote: >>>>> On Sun, 17 Mar 2013 19:31:38 -0400, Phil Hobbs >>>>> <pcdhSpamMeSenseless@electrooptical.net> wrote: >>>>> >>>>>> On 3/17/2013 6:44 PM, rickman wrote: >>>>>>> On 3/17/2013 6:20 PM, John Larkin wrote: >>>>>>>> >>>>>>>> It would be interesting to distort the table entries to minimize THD. The >>>>>>>> problem there is how to measure very low THD to close the loop. >>>>>>> >>>>>>> It could also be a difficult problem to solve without some sort of >>>>>>> exhaustive search. Each point you fudge affects two curve segments and >>>>>>> each curve segment is affected by two points. So there is some degree >>>>>>> of interaction. >>>>>>> >>>>>>> Anyone know if there is a method to solve such a problem easily? I am >>>>>>> sure that my current approach gets each segment end point to within &#2013266097;1 >>>>>>> lsb and I suppose once you measure THD in some way it would not be an >>>>>>> excessive amount of work to tune the 256 end points in a few passes >>>>>>> through. This sounds like a tree search with pruning. >>>>>>> >>>>>>> I'm assuming there would be no practical closed form solution. But >>>>>>> couldn't the THD be calculated in a simulation rather than having to be >>>>>>> measured on the bench? >>>>>>> >>>>>> >>>>>> 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. >>>>>> >>>>> >>>>> 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. >>>>> >>>>> >>>> >>>> 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. >>> >>> >> >>You're no fun anymore. :-) >> >>It does have a 0.01% FS range, so you'd expect to get a lot better than >> -80 dB. >> >>But for fast stuff, I'd probably do the same as you--filter the >>daylights out of it and amplify the residuals. >> >>Cheers >> >>Phil Hobbs > >Which is how most old school distortion analyzers work. Pop'tronics did a >four part series on one that used auto-tuned parallel-t filters with over >80 db fundamental notch rejection in the late 60s or early 70s. Realy a >slick design in a lot of ways, especially for a hobby mag.
That's how the HP339A works, with photoresistors used to fine tune the notches. There was an equivalent Heathkit.
> >?-)
It gets challenging at higher frequencies (MHz) and really low distortion levels. -- John Larkin Highland Technology Inc www.highlandtechnology.com jlarkin at highlandtechnology dot com Precision electronic instrumentation Picosecond-resolution Digital Delay and Pulse generators Custom timing and laser controllers Photonics and fiberoptic TTL data links VME analog, thermocouple, LVDT, synchro, tachometer Multichannel arbitrary waveform generators
On Wed, 20 Mar 2013 20:53:13 -0700, josephkk
<joseph_barrett@sbcglobal.net> wrote:

>On Sun, 17 Mar 2013 17:22:59 -0700, John Larkin ><jjlarkin@highNOTlandTHIStechnologyPART.com> wrote: > >> >>>I don't think you are clear on this. The table and linear interp are as >>>good as they can get with the exception of one or possibly two lsbs to >>>minimize the error in the linear interp. These is no way to correct for >>>this by adding a second harmonic, etc. Remember, this is not a pure >>>table lookup, it is a two step approach. >> >>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. > >How do you get that? All harmonics would be within the 90 degree table.
Second?
>Harmonic phase shifts is a different matter though, that would take >careful characterization of production quantities of the final product, >$$$$$.
Of course you'd have to null each harmonic for gain and phase, preferably on each production unit, but it wouldn't be difficult or expensive *if* you had a distortion analyzer with the required PPM performance. We regularly do automated calibrations, polynomial curve fits and such, more complex than this. -- John Larkin Highland Technology, Inc jlarkin at highlandtechnology dot com http://www.highlandtechnology.com Precision electronic instrumentation Picosecond-resolution Digital Delay and Pulse generators Custom laser drivers and controllers Photonics and fiberoptic TTL data links VME thermocouple, LVDT, synchro acquisition and simulation
Rickman:
To be clear, I was talking about the method to initialize the ROMs at synthesis-time, not at reset.  For this sort of thing, I really like the sfixed library (use ieee.fixed_pkg.all) and the math library (use ieee.math_real.all).  Then you have a function something like this:

  function fill_lut(lut_depth : integer;
                    nbits_out : integer)
    return sfixed_array is
    variable cos_lut : sfixed_vector(0 to lut_depth-1)(0 downto -(NBITS-1));
  begin
    for k in 0 to lut_depth-1 loop
      cos_lut(k) := 
        to_sfixed(cos(real(k)/real(lut_depth)*2.0*MATH_PI),cos_lut(0)); 
    end loop;
    return cos_lut;
  end function;

The sfixed_array is an array of sfixed, a type you have to declare yourself.  Then your actual ROM is a signal that is also an array of sfixed that you initialize by calling this function.  The values will be rounded using, I believe, a round-to-even convergent round.  MATH_PI is a constant that is built into the math library.
On 3/20/2013 8:11 PM, josephkk wrote:
> On Sun, 17 Mar 2013 17:56:21 -0400, rickman<gnuarm@gmail.com> wrote: > >> On 3/16/2013 8:02 PM, langwadt@fonz.dk wrote: >>> >>> why not skip the lut and just do a full sine approximation >>> >>> something like this, page 53-54 >>> >>> http://www.emt.uni-linz.ac.at/education/Inhalte/pr_mderf/downloads/ADSP-2181/Using%20the%20ADSP-2100%20Family%20Volume%201/Using%20the%20ADSP-2100%20Family%20Volume%201.pdf >> >> I'm not sure this would be easier. The LUT an interpolation require a >> reasonable amount of logic. This requires raising X to the 5th power >> and five constant multiplies. My FPGA doesn't have multipliers and this >> may be too much for the poor little chip. I suppose I could do some of >> this with a LUT and linear interpolation... lol > > No multipliers(?); ouch. My Simpson's rule method may fit but it will > need a trapezoidal strip out of a full Wallace tree multiplier. Preferably > three (speed issue, that may not be that big for you). How do you do the > linear interpolation without a multiplier?
I can do multiplies without a multiplier block. It gets to be problematic to do a *lot* of multiplies. The interp will be done in a shift/add type multiplier as well as the 8 bit amplitude multiplier for the scale factor (which is already in an existing product). The CODEC I am using has 128 clock cycle per sample so there is time to use a serial Booth's multiplier or two. -- Rick