Генерировать синусоидальный сигнал на C без использования стандартной функции

Ответ 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     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 для улучшения идей.

Ответ 2

... лучшая идея реализовать синусоидальный генератор в C?

Edit: предложите сначала прочитать эту статью, чтобы получить оценку того, что спрашивает OP.

Алгоритм CORDIC (Координированный поворот DIgital Computer) очень подходит для использования на небольших версиях 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). Следующая иллюстрация содержит представление фактических данных, полученных между голубыми точками. Остальная часть синусоидального сигнала создается посредством зеркалирования и переноса исходного раздела данных.

Ответ 3

Для создания точной синусоидальной функции требуется количество ресурсов (циклов ЦП и памяти), которые неоправданны в этом приложении. Ваша цель генерировать "гладкую" кривую синуса не учитывает требования приложения.

• Когда вы рисуете кривую, вы можете наблюдать несовершенства, когда вы применяете эту кривую к приводу LED PWM, человеческий глаз вообще не будет воспринимать эти недостатки.

• Также человеческий глаз не может воспринимать разницу в яркости между соседними значениями даже на 40-ступенчатой кривой, поэтому интерполяция не требуется.

• Это будет более эффективно, если вы создадите синусоидальную функцию, которая генерирует соответствующие значения PWM-диска напрямую без плавающей запятой. На самом деле, а не синусоидальной функции, более подходящий масштабированный поднятый косинус будет более подходящим, так что входной сигнал нуля приведет к выходу нуля, а вход половины количества значений в цикле приведет к максимальному значению для вашего PWM-диска,

Следующая функция генерирует кривую повышенного косинуса для 8-битного PWM FSD из 16 значений (и 16 байтов), генерирующих 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-разрядное разрешение более чем достаточно для задачи.

Ответ 4

Вы рассматривали моделирование части кривой синуса с [0..PI] в качестве параболы? Если яркость светодиода предназначена только для наблюдения человеческого глаза, формы кривых должны быть достаточно похожи, чтобы было обнаружено небольшое различие.

Вам просто нужно найти подходящее уравнение для его описания.

Хм...

Вершина при (PI/2, 1)

Пересечения по оси Х в (0, 0) и (PI, 0)

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

Где К будет...

K = 4 / (PI * PI)

Ответ 5

Классический хак, чтобы нарисовать круг (и, следовательно, генерировать синусоидальную волну тоже), - это 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, но вы можете настроить, разделив его на максимум или округление. Кроме того, можно использовать целочисленную арифметику, и вы можете исключить умножение/деление, используя отрицательную силу для двух для e, поэтому вместо этого можно использовать сдвиг.

Ответ 6

Для светодиода вы, вероятно, можете просто сделать с 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 с правильно определенной символической константой, обобщить функцию дальше и реализовать вычисление отрицательной фазы тоже).}

Ответ 7

Я бы пошел с приближением Бхаскары я синусоидальной функции. Используя градусы, от 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);
}

Вот график, показывающий точки, рассчитанные с этим приближением, а также точки, рассчитанные с помощью синусоидальной функции. Вы едва можете увидеть, как точки аппроксимации выглядывают из-под реальных точек синуса.

Ответ 8

Если ваше приложение не требует реальной точности, не убейте себя, придумав алгоритм 40-сильной или косинусоидальной волны. Кроме того, значения в вашей таблице должны соответствовать диапазону вашего светодиодного входа pwm.

Тем не менее, я взглянул на ваш код и вывел его и понял, что вы не интерполировали точки. С небольшим изменением я исправил его, и ошибка между функцией знака excel и вашим отключена на максимум около 0.0032 или около того. Это изменение довольно легко реализовать и было протестировано с использованием tcc, моего личного решения для тестирования алгоритма C.

Во-первых, я добавил еще один аргумент в ваш синусоидальный массив. Последняя точка имеет то же значение, что и первый элемент массива sine. Это фиксирует математику в вашей синусоидальной функции, в частности, когда вы устанавливаете x1 в (int) фазу% 40 и x2 на x1 + 1. Добавление дополнительной точки не требуется, поскольку вы можете установить x2 в (x1 + 1)% 40, но я выбрал первый подход. Я просто указываю на то, как вы могли бы это сделать. Я также добавил расчет остатка (в основном фаза - (int)). Я использую остаток для интерполяции. Я также добавил временный держатель стоимости синуса и переменную дельты.

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

Результаты выглядят неплохо. Сравнивая линейную аппроксимацию vs с системной функцией синуса с плавающей запятой, я дал график ошибок, показанный ниже.

Ответ 9

Вы могли бы использовать первые несколько терминов расширения серии Тейлора sin. Вы можете использовать как можно меньше терминов для достижения вашего предполагаемого уровня точности - несколько терминов, чем приведенный ниже пример, должны начинать сталкиваться с ограничениями 32-битного поплавка.

Пример:

#include <stdio.h>

// Please use the built-in floor function if you can. 
float my_floor(float f) {
    return (float) (int) f;
}

// Please use the built-in fmod function if you can.
float my_fmod(float f, float n) {
    return f - n * my_floor(f / n);
}

// t should be in given in radians.
float sin_t(float t) {
    const float PI = 3.14159265359f;

    // First we clamp t to the interval [0, 2*pi) 
    // because this approximation loses precision for 
    // values of t not close to 0. We do this by 
    // taking fmod(t, 2*pi) because sin is a periodic
    // function with period 2*pi.
    t = my_fmod(t, 2.0f * PI);

    // Next we clamp to [-pi, pi] to get our t as
    // close to 0 as possible. We "reflect" any values 
    // greater than pi by subtracting them from pi. This 
    // works because sin is an odd function and so 
    // sin(-t) = -sin(t), and the particular shape of sin
    // combined with the choice of pi as the endpoint
    // takes care of the negative.
    if (t >= PI) {
        t = PI - t;
    }

    // You can precompute these if you want, but
    // the compiler will probably optimize them out.
    // These are the reciprocals of odd factorials.
    // (1/n! for odd n)
    const float c0 = 1.0f;
    const float c1 = c0 / (2.0f * 3.0f);
    const float c2 = c1 / (4.0f * 5.0f);
    const float c3 = c2 / (6.0f * 7.0f);
    const float c4 = c3 / (8.0f * 9.0f);
    const float c5 = c4 / (10.0f * 11.0f);
    const float c6 = c5 / (12.0f * 13.0f);
    const float c7 = c6 / (14.0f * 15.0f);
    const float c8 = c7 / (16.0f * 17.0f);

    // Increasing odd powers of t.
    const float t3  = t * t * t;
    const float t5  = t3 * t * t;
    const float t7  = t5 * t * t;
    const float t9  = t7 * t * t;
    const float t11 = t9 * t * t;
    const float t13 = t9 * t * t;
    const float t15 = t9 * t * t;
    const float t17 = t9 * t * t;

    return c0 * t - c1 * t3 + c2 * t5 - c3 * t7 + c4 * t9 - c5 * t11 + c6 * t13 - c7 * t15 + c8 * t17;
}

// Test the output
int main() {
    const float PI = 3.14159265359f;
    float t;

    for (t = 0.0f; t < 12.0f * PI; t += (PI * 0.25f)) {
        printf("sin(%f) = %f\n", t, sin_t(t));
    }

    return 0;
}

Пример вывода:

sin(0.000000) = 0.000000
sin(0.785398) = 0.707107
sin(1.570796) = 1.000000
sin(2.356194) = 0.707098
sin(3.141593) = 0.000000
sin(3.926991) = -0.707107
sin(4.712389) = -1.000000
sin(5.497787) = -0.707098
sin(6.283185) = 0.000398
...
sin(31.415936) = 0.000008
sin(32.201332) = 0.707111
sin(32.986729) = 1.000000
sin(33.772125) = 0.707096
sin(34.557522) = -0.000001
sin(35.342918) = -0.707106
sin(36.128315) = -1.000000
sin(36.913712) = -0.707100
sin(37.699108) = 0.000393

Как вы можете видеть, есть еще кое-что для улучшения точности. Я не гений с арифметикой с плавающей запятой, поэтому, вероятно, некоторые из них связаны с реализацией floor/fmod или определенным порядком, в котором выполняются математические операции.

Ответ 10

Поскольку вы пытаетесь генерировать сигнал, я думаю, что использование дифференциального уравнения не должно быть плохой идеей! он дает что-то подобное

#include <stdlib.h>
#include <stdio.h>

#define DT (0.01f) //1/s
#define W0 (3)     //rad/s

int main(void) {
    float a = 0.0f;
    float b = DT * W0;
    float tmp;

    for (int i = 0; i < 400; i++) {
        tmp = (1 / (1 + (DT * DT * W0 * W0))) * (2 * a - b);
        b = a;
        a = tmp;

        printf("%f\n", tmp);
    }
}

Все еще установленная амплитуда и частота сигнала - это боль в шее:/

Ответ 11

Это поможет, если вы объясните, почему вы не хотите встроенной функции, но, как говорили другие, серия Taylor - это один из способов оценить значение. Тем не менее, другие ответы, похоже, фактически используют серию Маклорена, а не Тейлор. У вас должна быть таблица поиска как синуса, так и косинуса. Затем найдите x ₀, ближайшее значение x в вашей таблице поиска, к x, которое вы хотите, и найдите d = xx ₀. затем

Sin (х) = Sin (х ₀₎ +cos (х ₀₎ * D-син (х ₀₎ * d ^2/2-сов (х ₀₎ * d ^3/6 +...

Если ваша таблица поиска такова, что d <.01, то вы получите более двух цифр точности за семестр.

Другой метод - использовать тот факт, что если x = x ₀ +d, то

sin (x) = sin (x ₀) * cos (d) +cos (x ₀) * sin (d)

Вы можете использовать таблицу поиска, чтобы получить sin (x ₀) и cos (x ₀), а затем использовать числа Maclaurin для получения cos (d) и sin (d).