11

I am sampling a sine wave at 48 kHz, the frequency range of my sine wave can vary from 0 to 20000 Hz with a step of about 100 Hz. I am using a lookup table approach. So I generate 4096 samples for a sine wave for 4096 different phases. I think the general idea behind this to increment the step size and use different step sizes for different frequncy. So I do the following (pseudo code). But I am not sure how the step size is going to be related to the frequency I want to generate the samples of the sine wave of? For example if my frequency is 15000 Hz what would be the step size that I have to traverse? Is my sample size (4096) too low for this?

 // Pseudocode
 uint16_t audio_sample[4096] = {...};
 NSTEP = freq; //???How is the step size going to be related to the freq here

 for(int i = 0; i < 4096; i = i+NSTEP)
 {
     sine_f(i) = audio_sample[i];
 }

Thanks in advance.

5 Answers 5

17

You're on the right track - first we need to generate a sine wave LUT:

const int Fs = 48000;       // sample rate (Hz)
const int LUT_SIZE = 1024;  // lookup table size

int16_t LUT[LUT_SIZE];      // our sine wave LUT

for (int i = 0; i < LUT_SIZE; ++i)
{
    LUT[i] = (int16_t)roundf(SHRT_MAX * sinf(2.0f * M_PI * (float)i / LUT_SIZE));
}                           // fill LUT with 16 bit sine wave sample values

Note that we only need to generate this LUT once, e.g. during initialisation.

Now that we have a sine wave LUT we can use it to generate any frequency we wish to using a phase accumulator:

const int BUFF_SIZE = 4096;  // size of output buffer (samples)
int16_t buff[BUFF_SIZE];     // output buffer

const int f = 1000;          // frequency we want to generate (Hz)

const float delta_phi = (float) f / Fs * LUT_SIZE;
                             // phase increment

float phase = 0.0f;          // phase accumulator

// generate buffer of output
for (int i = 0; i < BUFF_SIZE; ++i)
{
    int phase_i = (int)phase;        // get integer part of our phase
    buff[i] = LUT[phase_i];          // get sample value from LUT
    phase += delta_phi;              // increment phase
    if (phase >= (float)LUT_SIZE)    // handle wraparound
        phase -= (float)LUT_SIZE;
}

Note: for higher quality output you can use linear interpolation between the LUT values at phase_i and phase_i + 1, but the above approach is good enough for most audio applications.

5
  • Thanks a lot for your anser. I got most parts except the way you are handling wraparound, when you do phase -=(float)LUT_SIZE I am confused it takes the phase back to the (3*pi/2,pi) interval not (0,pi/2) as I would imagine. Can you please explain the wraparound part a bit? Again thanks so much! Nov 20, 2012 at 6:11
  • 2
    I think you're getting confused between the LUT and the output buffer - I deliberately made them different sizes to illustrate the fact that they are independent and unrelated. When the phase exceeds the size of the LUT we make it wrap around by subtracting LUT_SIZE.
    – Paul R
    Nov 20, 2012 at 6:14
  • Thanks! You are right I confused them. so that means samples 1024-4096 will all be mapped to sin(2*pi), right? Nov 20, 2012 at 6:24
  • No - it depends on the frequency - we are stepping through the LUT with an increment of delta_phi - lower frequencies have a smaller delta_phi, higher frequencies have a higher delta_phi - try plugging in some values to see how this works.
    – Paul R
    Nov 20, 2012 at 6:27
  • 1
    +1 for answering the question well AND mentioning linear interpolation! Nov 20, 2012 at 14:03
5

With look up table approach one may wish to use the memory efficiently and store only the first quadrant of the sine wave.

Then second quadrant = sin[180 - angle]  ; // for angles 90..180 degrees
third quadrant = -sin[angle-180];          // for 180..270
fourth quadrant = -sin[-angle+360]         // for 270..360

I would recommend linear interpolation, but there's also vector rotation based approach for sine generation (that produces both sine and cosine simultaneously)

x_0 = 1, y_0 = 0;
x_n+1 = x_n * cos(alpha) - y_n * sin(alpha)
y_n+1 = x_n * sin(alpha) + y_n * cos(alpha), 

where alpha=phase difference of the frequency == 2pi*fHz/Ts, with fHz is frequency to be produced (in Hertz) and Ts is sampling time (or 1/Ts = sampling frequenzy eg. 44100 Hz).

which leads to a resonating feedback filter approach, whose transfer function f(z) has two conjugate poles at unit circle (z=e^jomegaT).

y[n] = y[n-1] - 2*cos(alpha) * y[n-2], with
y[0] = 0,
y[1] = sin(alpha)

The fun part is that one can change the alpha (cos(alpha)) on the fly. The downside of this IIR filter approach is that it's unstable by definition. Floating point and especially fixed point inaccuracies accumulate and lead to either exponential decay, or exponential growth of the magnitude. That can however be cured with allowing a slight phase distortion.

Instead of as in CORDIC rotation having a known per iteration amplification factor:

K = sqrt(1+sin(alpha)^2);

x' = x - y*sin(alpha);
y' = y + x*sin(alpha);

one can do

x' = x - y * sin(alpha);
y'' = y + x' * sin(alpha);

which doesn't produce perfect circles for (x', y'') but stable ellipses even with fixed point arithmetic. (Note that this assumes relatively small values of alpha, meaning also relatively low frequencies.)

3

If you want high precision, you can use a trig identities to have both a small LUT, and clean sine waves.

static float gTrigSinHi[256], gTrigSinLo[256];
static float gTrigCosHi[256], gTrigCosLo[256];

////////////////////////////////////////
// Sets up the fast trig functions
void FastTrigInit()
{
    unsigned i;
    for(i = 0; i < 256; ++i)
    {
        gTrigSinHi[i] = sin(2.0 * M_PI / 0xFFFF * (i << 8));
        gTrigSinLo[i] = sin(2.0 * M_PI / 0xFFFF * (i << 0));
        gTrigCosHi[i] = cos(2.0 * M_PI / 0xFFFF * (i << 8));
        gTrigCosLo[i] = cos(2.0 * M_PI / 0xFFFF * (i << 0));
    }
}

////////////////////////////////////////
// Implements sin as
//      sin(u+v) = sin(u)*cos(v) + cos(u)*sin(v)
float FastSin(unsigned short val)
{
    unsigned char hi = (val >> 8);
    unsigned char lo = (val & 0xFF);
    return gTrigSinHi[hi] * gTrigCosLo[lo]
        + gTrigCosHi[hi] * gTrigSinLo[lo];
}

////////////////////////////////////////
// Implements cos as
//      cos(u+v) = cos(u)*cos(v) - sin(u)*sin(v)
float FastCos(unsigned short val)
{
    unsigned char hi = (val >> 8);
    unsigned char lo = (val & 0xFF);
    return gTrigCosHi[hi] * gTrigCosLo[lo]
        - gTrigSinHi[hi] * gTrigSinLo[lo];
}
1
  • 2
    That's a good one. And it can be reduced with cos(n)=sin[(n+64)&255]; plus the phase accumulation is easy to do with integers (or 16.16 fixed points) Nov 20, 2012 at 18:18
0

Very good answer, this is classic software DDS. Facing the same problem these days. There is no need to use floats

        UInt16 f = 400; // Hz, up to 16384 :)
        UInt16 delta = (UInt16)(((UInt32)f * LUT_SIZE ) / fmt_SamplesPerSec);
        UInt16 phase = 0;

        for (int i = 0; i < BUFF_SIZE; ++i)
        {
            buff[i] = LUT[phase];
            phase += delta;
            phase &= LUT_SIZE-1;
        }

Let phase wraparound LUT size as mask. And don't care about using quadrants since for my purpose I have a huge MIPS for this requirements already.

-1

According to "http://en.wikipedia.org/wiki/X86_instruction_listings" if you have x80387 or newer there is a sine instruction, so just call it directly. You just have to figure out how to add some in-line assembly language to your program. This way you don't have to worry if your input value is not an exact match to what is in your table.

1
  • 1
    I think if the poster didn't want to use a lookup table, the poster would use sinf() before resorting to assembly. Especially since you can coerce sinf() to compile to the sine instruction on many compilers. Nov 20, 2012 at 5:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.