How to look up sine of different frequencies from

问题:

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.

回答1:

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.

回答2:

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];
}

回答4:

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.

回答5:

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.