如何从固定大小的查找表中查找不同频率的正弦波？

Question

我正在对 48 kHz 的正弦波进行采样，我的正弦波的频率范围可以从 0 到 20000 Hz 变化，步长约为 100 Hz。我正在使用查找表方法。因此，我为 4096 个不同相位的正弦波生成 4096 个样本。我认为这背后的一般想法是增加步长并为不同的频率使用不同的步长。所以我做了以下（伪代码）。但我不确定步长与我想要生成正弦波样本的频率有什么关系？例如，如果我的频率是 15000 Hz，那么我必须遍历的步长是多少？我的样本量（4096）是否太低了？

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

提前致谢。

Answer 1

您在正确的轨道上 - 首先我们需要生成一个正弦波 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

请注意，我们只需要生成此 LUT 一次，例如在初始化期间。

现在我们有了一个正弦波 LUT，我们可以使用它来生成我们希望使用 相位累加器的任何频率：

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

注意：为了获得更高质量的输出，您可以在 phase_i 和 phase_i + 1 的 LUT 值之间使用线性插值，但上述方法对于大多数音频应用来说已经足够了。

Answer 2

很好的答案，这是经典软件DDS。 这些天面临同样的问题。不需要使用浮点数

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

让相位环绕 LUT 大小作为掩码。并且不要在意使用象限，因为出于我的目的，我已经有一个巨大的 MIPS 来满足这个要求。

Answer 3

如果您想要高精度，您可以使用三角标识来同时拥有小的 LUT 和干净的正弦波。

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

Answer 4

使用查找表方法，您可能希望有效地使用内存并仅存储正弦波的第一象限。

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

我会推荐线性插值，但也有基于矢量旋转的正弦生成方法（同时产生正弦和余弦）

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), 

其中 alpha=频率的相位差 == 2pi*fHz/Ts，其中 fHz 是要产生的频率（以赫兹为单位），Ts 是采样时间（或 1/Ts = 采样频率例如 44100 赫兹）。

这导致了一种谐振反馈滤波器方法，其传递函数 f(z) 在单位圆 (z=e^jomegaT) 处有两个共轭极点。

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

有趣的是，人们可以随时更改 alpha (cos(alpha))。这种 IIR 滤波器方法的缺点是根据定义它是不稳定的。浮点，尤其是定点的不准确性累积并导致指数衰减或数量级的指数增长。然而，这可以通过允许轻微的相位失真来解决。

而不是 CORDIC 旋转具有已知的每次迭代放大系数：

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);

它不会为 (x', y'') 生成完美的圆，但即使使用定点算术也可以生成稳定的椭圆。 （请注意，这假设 alpha 值相对较小，也意味着频率相对较低。）

Answer 5

根据“http://en.wikipedia.org/wiki/X86_instruction_listings”如果你有 x80387 或更新版本有一个正弦指令，所以直接调用它。您只需要弄清楚如何将一些内联汇编语言添加到您的程序中。这样您就不必担心输入值与表中的值不完全匹配。