在不使用标准功能的情况下在C中生成正弦信号

Question

我想在不使用标准函数sin()的情况下在C中生成正弦信号,以触发LED亮度的正弦形状变化.我的基本想法是使用具有40个点和插值的查找表.

这是我的第一个方法:

const int sine_table[40] = {0, 5125, 10125, 14876, 19260, 23170, 26509, 29196,
31163, 32364, 32767,  32364, 31163, 29196, 26509, 23170, 19260, 14876, 10125,
5125, 0, -5126, -10126,-14877, -19261, -23171, -26510, -29197, -31164, -32365,
-32768, -32365, -31164, -29197, -26510, -23171, -19261, -14877, -10126, -5126};

int i = 0;
int x1 = 0;
int x2 = 0;
float y = 0;

float sin1(float phase)
{
    x1 = (int) phase % 41;
    x2 = x1 + 1;
    y = (sine_table[x2] - sine_table[x1])*((float) ((int) (40*0.001*i*100) % 4100)/100 - x1) + sine_table[x1];
    return y;
}

int main()
{
    while(1)
    {
    printf("%f      ", sin1(40*0.001*i)/32768);
    i = i + 1;
    }
}

不幸的是,这个功能有时会返回值远小于1大.此外,插值似乎并不好(我用它来创建一个LED的正弦波形的亮度变化,但这些都是非常unsmoooth).

有人有更好的想法在C中实现正弦发生器吗？

Answer 1

OP的主要问题是为表查找生成索引.

OP的代码尝试访问外部数组sine_table[40]导致未定义的行为.至少修复一下.

const int sine_table[40] = {0, 5125, 10125, ...
    ...
    x1 = (int) phase % 41;                     // -40 <= x1 <= 40
    x2 = x1 + 1;                               // -39 <= x2 <= 41  
    y = (sine_table[x2] - sine_table[x1])*...  // bad code, consider x1 = 40 or x2 = 40,41

建议的改变

    x1 = (int) phase % 40;   // mod 40, not 41
    if (x1 < 0) x1 += 40;    // Handle negative values
    x2 = (x1 + 1) % 40;      // Handle wrap-around 
    y = (sine_table[x2] - sine_table[x1])*...  

有更好的方法,但关注OP的方法,见下文.

#include <math.h>
#include <stdio.h>

const int sine_table[40] = { 0, 5125, 10125, 14876, 19260, 23170, 26509, 29196,
31163, 32364, 32767, 32364, 31163, 29196, 26509, 23170, 19260, 14876, 10125,
5125, 0, -5126, -10126, -14877, -19261, -23171, -26510, -29197, -31164, -32365,
-32768, -32365, -31164, -29197, -26510, -23171, -19261, -14877, -10126, -5126 };

int i = 0;
int x1 = 0;
int x2 = 0;
float y = 0;

float sin1(float phase) {
  x1 = (int) phase % 40;
  if (x1 < 0) x1 += 40;
  x2 = (x1 + 1) % 40;
  y = (sine_table[x2] - sine_table[x1])
      * ((float) ((int) (40 * 0.001 * i * 100) % 4100) / 100 - x1)
      + sine_table[x1];
  return y;
}

int main(void) {
  double pi = 3.1415926535897932384626433832795;
  for (int j = 0; j < 1000; j++) {
    float x = 40 * 0.001 * i;
    float radians = x * 2 * pi / 40;
    printf("%f %f %f\n", x, sin1(x) / 32768, sin(radians));
    i = i + 1;
  }
}

产量

         OP's     Reference sin()
0.000000 0.000000 0.000000
0.040000 0.006256 0.006283
0.080000 0.012512 0.012566
...
1.960000 0.301361 0.303035
2.000000 0.308990 0.309017
2.040000 0.314790 0.314987
...
39.880001 -0.020336 -0.018848
39.919998 -0.014079 -0.012567
39.959999 -0.006257 -0.006283

更好的代码不会将值i, x1, x2, y作为全局变量传递,而是作为函数参数或函数变量传递.也许这是OP调试的神器.

有人有更好的想法在C中实现正弦发生器吗？

这是相当广泛的.速度,精度,代码空间,可移植性或可维护性更好？ sine()功能很容易.高质量的人需要付出更多努力.

虽然模糊,OP使用小型查找表是一个很好的开始 - 虽然我看到它可以在没有任何浮点数学的情况下完成.我建议OP构建一个经过测试和运行的解决方案,并将其发布在Code Review中以获得改进的想法.

Answer 2

...在C中实现正弦发生器更好的主意？

编辑:建议先阅读本文,以了解OP的要求.

从您问题中提供的上下文中,我假设更好的单词可能与编译代码的大小和速度有关,可能需要在小型微处理器上运行.

CORDIC(坐标旋转数字计算机)算法非常适用于较小的uP和具有有限数学计算能力的FPGA实现,因为它仅使用基本算术(加法,减法和移位)计算值的正弦和余弦.更多关于CORDIC,以及如何使用它来产生一个角度的正弦/余弦这里提供.

还有几个站点提供算法实现示例. 简单CORDIC包含有关如何生成表格的详细说明,该表格可以预编译以便在目标设备上使用,以及用于测试以下函数输出的代码(使用定点数学运算):

(参见以下文档和链接中的其他功能)

#define cordic_1K 0x26DD3B6A
#define half_pi 0x6487ED51
#define MUL 1073741824.000000
#define CORDIC_NTAB 32
int cordic_ctab [] = {0x3243F6A8, 0x1DAC6705, 0x0FADBAFC, 0x07F56EA6, 0x03FEAB76, 0x01FFD55B, 
0x00FFFAAA, 0x007FFF55, 0x003FFFEA, 0x001FFFFD, 0x000FFFFF, 0x0007FFFF, 0x0003FFFF, 
0x0001FFFF, 0x0000FFFF, 0x00007FFF, 0x00003FFF, 0x00001FFF, 0x00000FFF, 0x000007FF, 
0x000003FF, 0x000001FF, 0x000000FF, 0x0000007F, 0x0000003F, 0x0000001F, 0x0000000F, 
0x00000008, 0x00000004, 0x00000002, 0x00000001, 0x00000000 };

void cordic(int theta, int *s, int *c, int n)
{
  int k, d, tx, ty, tz;
  int x=cordic_1K,y=0,z=theta;
  n = (n>CORDIC_NTAB) ? CORDIC_NTAB : n;
  for (k=0; k<n; ++k)
  {
    d = z>>31;
    //get sign. for other architectures, you might want to use the more portable version
    //d = z>=0 ? 0 : -1;
    tx = x - (((y>>k) ^ d) - d);
    ty = y + (((x>>k) ^ d) - d);
    tz = z - ((cordic_ctab[k] ^ d) - d);
    x = tx; y = ty; z = tz;
  }  
 *c = x; *s = y;
}

编辑:
我发现使用 Simple CORDIC站点上的示例的文档非常容易理解.但是,我遇到的一件小事是在编译文件cordic-test.c错误时:使用未声明的标识符'M_PI'.看来,当执行编译gentable.c文件(生成cordic-test.c文件)时,行:

#define M_PI 3.1415926535897932384626

虽然包含在自己的声明中,但未包含在用于生成文件的printf语句中cordic-test.c.一旦解决了这个问题,一切都按照宣传的方式进行.

如记载的那样,产生的数据范围产生1/4的完整正弦周期(-π/ 2-π/ 2).下图包含浅蓝点之间产生的实际数据的表示.通过镜像和转置原始数据部分来制造正弦信号的其余部分.

Answer 3

生成准确的正弦函数需要一定量的资源(CPU周期和内存),这在此应用程序中是没有根据的.您生成"平滑"正弦曲线的目的是未能考虑应用的要求.

• 当您绘制曲线时,您可能会观察到瑕疵,当您将该曲线应用于LED PWM驱动器时,人眼根本不会察觉到那些瑕疵.

• 即使是40步曲线,人眼也不可能察觉到相邻值之间的亮度差异,因此不需要插值.

• 如果生成正弦函数直接生成没有浮点的适当PWM驱动值,则效率会更高.事实上,不是正弦函数,缩放的升余弦更合适,因此输入为零会导致输出为零,并且循环中值的一半输入会导致PWM驱动的最大值.

以下函数从16位(和16字节)查找生成8位FSD PWM的升余弦曲线,生成59步循环.因此,与40步浮点实现相比,它具有内存和性能效率.

#include <stdint.h>

#define LOOKUP_SIZE 16
#define PWM_THETA_MAX (LOOKUP_SIZE * 4 - 4)

uint8_t RaisedCosine8bit( unsigned n )
{
    static const uint8_t lookup[LOOKUP_SIZE] = { 0, 1, 5, 9,
                                                 14, 21, 28, 36,
                                                 46, 56, 67, 78,
                                                 90, 102, 114, 127} ;
    uint8_t s = 0 ;
    n = n % PWM_THETA_MAX ;

    if( n < LOOKUP_SIZE )
    {
        s = lookup[n] ;
    }
    else if( n < LOOKUP_SIZE * 2 - 1 )
    {
        s = 255 - lookup[LOOKUP_SIZE * 2 - n - 2] ;
    }
    else if( n < LOOKUP_SIZE * 3 - 2 )
    {
        s = 255 - lookup[n - LOOKUP_SIZE * 2 + 2] ;
    }
    else
    {
        s = lookup[LOOKUP_SIZE * 4 - n - 4] ;
    }

    return s ;
}

对于0 <= theta < PWM_THETA_MAX曲线的输入,如下所示:

我建议足够光滑照明.

在实践中,您可以使用它:

for(;;)
{
    for( unsigned i = 0; i < PWM_THETA_MAX; i++ )
    {
        LedPwmDrive( RaisedCosine8bit( i ) ) ;
        Delay( LED_UPDATE_DLEAY ) ;
    }
}

如果PWM范围不是0到255,只需缩放功能的输出; 8位分辨率足以完成任务.

Answer 4

绘制圆圈的经典黑客(因此也产生正弦波)是Marvin Minsky的Hakmem#149.例如,:

#include <stdio.h>

int main(void)
{
    float x = 1, y = 0;

    const float e = .04;

    for (int i = 0; i < 100; ++i)
    {
        x -= e*y;
        y += e*x;
        printf("%g\n", y);
    }
}

它会略微偏心,不是一个完美的圆圈,你可能会略微超过1,但你可以通过除以最大值或舍入来调整.此外,可以使用整数运算,并且可以通过使用2的负幂来消除乘法/除法e,因此可以使用shift.

Answer 5

您是否考虑将[0..PI]的正弦曲线部分建模为抛物线？如果LED的亮度仅用于人眼观察,则曲线的形状应该足够相似,以便检测到很小的差异.

您只需要找出适当的等式来描述它.

嗯,......

顶点(PI/2,1)

(0,0)和(PI,0)处的X轴交点

f(x) = 1 - K * (x - PI/2) * (x - PI/2)

K会是......

K = 4 / (PI * PI)

Answer 6

对于LED,你可能只需要16步左右而不进行插值.也就是说,我可以在你的sin1()函数中看到至少两个奇怪的东西:

1)您有40个数据点sine_table,但您正在x1使用输入的索引模41.这似乎不是处理周期性的正确方法,并且让我们x1指向数组的最后一个索引.
2)然后你添加+1,所以x2可以超过数组的限制.
3)您正在使用i该功能,但它仅在主程序中设置.我不知道它应该做什么,但在一个简单的计算函数中使用这样的全局似乎很脏.也许它应该为插值提供小数部分,但不应该使用phase它.

这是一个简单的插值器,似乎有效.调整味道.

#include <assert.h>

int A[4] = {100, 200, 400, 800};    
int interpolate(float x)
{
    if (x == 3.00) {
        return A[3];
    }
    if (x > 3) {
        return interpolate(6 - x);
    }
    assert(x >= 0 && x < 3);
    int i = x;
    float frac = x - i;
    return A[i] + frac * (A[i+1] - A[i]);
}

一些任意的样本输出:

interpolate(0.000000) = 100
interpolate(0.250000) = 125
interpolate(0.500000) = 150
interpolate(1.000000) = 200
interpolate(1.500000) = 300
interpolate(2.250000) = 500
interpolate(2.999900) = 799
interpolate(3.000000) = 800
interpolate(3.750000) = 500

_{(我将把它留给感兴趣的读者3用一个正确定义的符号常量替换所有出现的,进一步概括函数,并实现计算负相位.)}

Answer 7

除非你的应用程序要求真正的精确度,否则不要自己想出一个40点正弦波或余弦波算法.此外,表格中的值应与LED的pwm输入范围相匹配.

也就是说,我看了你的代码和它的输出,并认为你没有插入点之间.经过一点修改,我修复了它,excel的sign函数和你之间的错误最多约为0.0032左右.这个变化很容易实现,并且已经使用tcc进行了测试,tcc是我个人的C算法测试的首选.

首先,我向你的正弦数组添加了一个点.最后一个点设置为与正弦数组中的第一个元素相同的值.这将修复正弦函数中的数学运算,特别是当您将x1设置为(int)phase%40,x2设置为x1 + 1时.添加额外点不是必需的,因为您可以将x2设置为(x1 + 1)%40,但我选择了第一种方法.我只想指出你可以通过不同的方式实现这一目标.我还添加了余数的计算(基本相 - (int)阶段).我正在使用余数进行插值.我还添加了一个临时正弦值持有者和一个delta变量.

const int sine_table[41] = 
{0, 5125, 10125, 14876, 19260, 23170, 26509, 29196,
31163, 32364, 32767,  32364, 31163, 29196, 26509, 23170, 
19260, 14876, 10125, 5125, 0, -5126, -10126,-14877,
-19261, -23171, -26510, -29197, -31164, -32365, -32768, -32365,
-31164, -29197, -26510, -23171, -19261, -14877, -10126, -5126, 0};

int i = 0;
int x1 = 0;
int x2 = 0;
float y = 0;

float sin1(float phase)
{
    int tsv,delta;
    float rem;

    rem = phase - (int)phase;
    x1 = (int) phase % 40;
    x2 = (x1 + 1);

    tsv=sine_table[x1];
    delta=sine_table[x2]-tsv;

    y = tsv + (int)(rem*delta);
    return y;
}

int main()
{
    int i;  
    for(i=0;i<420;i++)
    {
       printf("%.2f, %f\n",0.1*i,sin1(0.1*i)/32768);
    }
    return 0;
}

结果看起来很不错.比较线性近似与系统的浮点正弦函数给出了如下所示的误差图.

Answer 8

我会选择Bhaskara,我接近正弦函数.使用度数,从0到180,您可以像这样近似值

float Sine0to180(float phase)
{
    return (4.0f * phase) * (180.0f - phase) / (40500.0f - phase * (180.0f - phase));
}

如果你想考虑任何角度,你会添加

float sine(float phase)
{
    float FactorFor180to360 = -1 * (((int) phase / 180) % 2 );
    float AbsoluteSineValue = Sine0to180(phase - (float)(180 * (int)(phase/180)));
    return AbsoluteSineValue * FactorFor180to360;
}

如果你想用弧度来做,你会添加

float SineRads(float phase)
{
    return Sine(phase * 180.0f / 3.1416);
}

下面的图表显示了使用此近似值计算的点数以及使用正弦函数计算的点数.你几乎看不到从实际正弦点下面偷看的近似点.