Реализация тригонометрических функций для приложений реального времени

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

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

Справочная таблица

#define TABLE_SIZE 256
#define PI         3.14
#define STEP_SIZE  (PI/(2*(TABLE_SIZE-1))) // N points divide PI/2 into N-1 segments

// look-up table containing values of sine function over one quarter of period
// angel step size is pi/(2*256) i.e. pi/512 approximately 0.35°
double sinLUT[TABLE_SIZE] = 
{
0.000000,  0.006160,  0.012320,  0.018479,  0.024637,  0.030795,  0.036951,  0.043107,  
0.049260,  0.055411,  0.061561,  0.067708,  0.073853,  0.079994,  0.086133,  0.092268,  
0.098400,  0.104528,  0.110653,  0.116773,  0.122888,  0.128999,  0.135105,  0.141206,  
0.147302,  0.153392,  0.159476,  0.165554,  0.171626,  0.177691,  0.183750,  0.189801,  
0.195845,  0.201882,  0.207912,  0.213933,  0.219946,  0.225951,  0.231948,  0.237935,  
0.243914,  0.249883,  0.255843,  0.261793,  0.267733,  0.273663,  0.279583,  0.285492,  
0.291390,  0.297277,  0.303153,  0.309017,  0.314870,  0.320710,  0.326539,  0.332355,  
0.338158,  0.343949,  0.349727,  0.355491,  0.361242,  0.366979,  0.372702,  0.378411,  
0.384106,  0.389786,  0.395451,  0.401102,  0.406737,  0.412356,  0.417960,  0.423549,  
0.429121,  0.434676,  0.440216,  0.445738,  0.451244,  0.456733,  0.462204,  0.467658,  
0.473094,  0.478512,  0.483911,  0.489293,  0.494656,  0.500000,  0.505325,  0.510631,  
0.515918,  0.521185,  0.526432,  0.531659,  0.536867,  0.542053,  0.547220,  0.552365,  
0.557489,  0.562593,  0.567675,  0.572735,  0.577774,  0.582791,  0.587785,  0.592758,  
0.597707,  0.602635,  0.607539,  0.612420,  0.617278,  0.622113,  0.626924,  0.631711,  
0.636474,  0.641213,  0.645928,  0.650618,  0.655284,  0.659925,  0.664540,  0.669131,  
0.673696,  0.678235,  0.682749,  0.687237,  0.691698,  0.696134,  0.700543,  0.704926,  
0.709281,  0.713610,  0.717912,  0.722186,  0.726434,  0.730653,  0.734845,  0.739009,  
0.743145,  0.747253,  0.751332,  0.755383,  0.759405,  0.763398,  0.767363,  0.771298,  
0.775204,  0.779081,  0.782928,  0.786745,  0.790532,  0.794290,  0.798017,  0.801714,  
0.805381,  0.809017,  0.812622,  0.816197,  0.819740,  0.823253,  0.826734,  0.830184,  
0.833602,  0.836989,  0.840344,  0.843667,  0.846958,  0.850217,  0.853444,  0.856638,  
0.859800,  0.862929,  0.866025,  0.869089,  0.872120,  0.875117,  0.878081,  0.881012,  
0.883910,  0.886774,  0.889604,  0.892401,  0.895163,  0.897892,  0.900587,  0.903247,  
0.905873,  0.908465,  0.911023,  0.913545,  0.916034,  0.918487,  0.920906,  0.923289,  
0.925638,  0.927951,  0.930229,  0.932472,  0.934680,  0.936852,  0.938988,  0.941089,  
0.943154,  0.945184,  0.947177,  0.949135,  0.951057,  0.952942,  0.954791,  0.956604,  
0.958381,  0.960122,  0.961826,  0.963493,  0.965124,  0.966718,  0.968276,  0.969797,  
0.971281,  0.972728,  0.974139,  0.975512,  0.976848,  0.978148,  0.979410,  0.980635,  
0.981823,  0.982973,  0.984086,  0.985162,  0.986201,  0.987202,  0.988165,  0.989092,  
0.989980,  0.990831,  0.991645,  0.992421,  0.993159,  0.993859,  0.994522,  0.995147,  
0.995734,  0.996284,  0.996795,  0.997269,  0.997705,  0.998103,  0.998464,  0.998786,  
0.999070,  0.999317,  0.999526,  0.999696,  0.999829,  0.999924,  0.999981,  1.000000
};

Расчет синусоидальной функции

double sine(double arg, double lut[TABLE_SIZE]){
double retval;
double rem;
uint8_t index;

if(arg >= 0 && arg <= PI/2){
    // first quadrant
    index = arg/STEP_SIZE;
    rem = arg - index*STEP_SIZE;
    if(rem > 0){
        // sine value for given argument isn't directly in the lut
        if(index == (TABLE_SIZE-1)){
            // last point in the lut so the interval for the interpolation
            // is the last interval bounded with the index-1 and index
            index -= 1;
        }
        retval = (lut[index+1] - lut[index])/STEP_SIZE*rem + lut[index];
    }else{
        // sine value for given argument is directly in the lut
        retval = lut[index];
    }
}else if(arg > PI/2 && arg <= PI){
    // second quadrant
    index = (PI - arg)/STEP_SIZE;
    rem = (PI - arg) - index*STEP_SIZE;
    if(rem > 0){
        if(index == (TABLE_SIZE-1)){
            index -= 1;
        }
        retval = (lut[index+1] - lut[index])/STEP_SIZE*rem + lut[index];
    }else{
        retval = lut[index];
    }
    
}else if(arg > PI && arg <= 3*PI/2){
    // third quadrant
    index = (arg - PI)/STEP_SIZE;
    rem = (arg - PI) - index*STEP_SIZE;
    if(rem > 0){
        if(index == (TABLE_SIZE-1)){
            index -= 1;
        }
        retval = (-lut[index+1] + lut[index])/STEP_SIZE*rem - lut[index];
    }else{
        retval = -lut[index];
    }
    
}else{
    // fourth quadrant
    index = (2*PI - arg)/STEP_SIZE;
    rem = (2*PI - arg) - index*STEP_SIZE;
    if(rem > 0){
        if(index == (TABLE_SIZE-1)){
            index -= 1;
        }
        retval = (-lut[index+1] + lut[index])/STEP_SIZE*rem - lut[index];
    }else{
        retval = -lut[index];
    }
}
return retval;
}

Функция косинуса

double cosine(double arg, double lut[TABLE_SIZE]){
 double temp;

 temp = (arg + PI/2);
 if(temp > 2*PI){
    temp -= 2*PI;
 }
 return sine(temp, lut);
}

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

Поэтому я обращаюсь не к тому, что вы сделали неправильно в своем коде, а к тому, что вы сделали неправильно, даже думая о своей реализации. Во-первых, погуглите, как реализовать косинус на микроконтроллере или как аппроксимировать косинус в целом. Абсолютно не попробуйте погуглить «Как реализовать X с помощью метода Y», потому что вы Только собираюсь получить этот метод (особенно не сначала). Это расширит то, что вы считали возможным, например, одним гораздо более быстрым способом, который много проще реализовать — использовать Чебышевское приближение. Вы бы заранее сгенерировали многочлены и скомпилировать их в в зависимости от того, насколько точным вы хотите, чтобы ваше приложение было. Вы также можете воспользоваться симметрией, чтобы рассчитывать только на каждые 45 градусов. Приближение Чебышева требует фиксированного количества времени, которое легко измерить.

Кроме того, вы можете не захотеть использовать PI для внутренних целей. Использовать базовый пи или же 2,0 * пи. Я имею в виду, когда вы генерируете многочлены Чебышева, делайте это так, чтобы ваш ввод был нормализованный 1.0 -> -1.0 или другой нормализованный диапазон. Это может помочь с проблемами точности (поскольку число Пи иррационально и в любом случае не может быть полностью представлено), а также позволяет использовать такие функции, как cos2pi и возвращать более точные результаты, если пользователь может легко определить угол в терминах нормализованного угла.

Другой метод, который намного быстрее, на самом деле использует реальное численное приближение, специфичное для функции аля https://stackoverflow.com/a/28050328/. Нет веток видишь? И есть дальнейшие усовершенствования подобных методов, которые могут дать еще более точные ответы за счет скорости. В целом эти методы превосходят чебышевскую формулировку как по скорости, так и по точности, но чебышевская формулировка не требует аналитического или же Для работы расширенной числовой формулировки для конкретного приложения нужны только входные и выходные значения, поэтому, если у вас есть более сложная функция, которая вас интересует только для определенного диапазона значений, ее часто можно использовать для аппроксимации всего этого.

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

• #define PI 3.14… можно было бы использовать еще несколько цифр!
• синус должен работать для чисел больше, чем 2.0 * PI.
• синус должен работать для отрицательных чисел.
• то же самое для косинуса.
• if(temp > 2*PI) { temp -= 2*PI; } неэффективен для чисел больше, чем 4.0 * PI.

if(rem > 0){
    // sine value for given argument isn't directly in the lut
    if(index == (TABLE_SIZE-1)){
        // last point in the lut so the interval for the interpolation
        // is the last interval bounded with the index-1 and index
        index -= 1;
    }
    retval = (lut[index+1] - lut[index])/STEP_SIZE*rem + lut[index];
}else{
    // sine value for given argument is directly in the lut
    retval = lut[index];
}

Действительно ли возможно получить остаток, если index == TABLE_SIZE-1?

Может, лучше было бы всегда сделать интерполяцию? Мы могли бы изменить таблицу на double sinLUT[TABLE_SIZE + 1] = ..., и добавьте дополнительный 1.0 значение до конца. Тогда у нас всегда есть index+1 доступно (примечание: нам придется изменить тип индекса с uint8_t).

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

Я был бы склонен сделать что-то вроде этого:

#include <assert.h>
#include <math.h>
#include <stdbool.h>
#include <stdint.h>

#define TABLE_SIZE 256
#define PI         3.1415926535897932384626433832795
#define STEP_SIZE  (1.0 / (TABLE_SIZE - 1))

// look-up table containing values of sine function over one quarter of period
// angle step size is pi/(2*256) i.e. pi/512 approximately 0.35°
double sinLUT[TABLE_SIZE + 1] = 
{
    0.000000,  0.006160,  0.012320,  0.018479,  0.024637,  0.030795,  0.036951,  0.043107, 
    0.049260,  0.055411,  0.061561,  0.067708,  0.073853,  0.079994,  0.086133,  0.092268, 
    0.098400,  0.104528,  0.110653,  0.116773,  0.122888,  0.128999,  0.135105,  0.141206, 
    0.147302,  0.153392,  0.159476,  0.165554,  0.171626,  0.177691,  0.183750,  0.189801, 
    0.195845,  0.201882,  0.207912,  0.213933,  0.219946,  0.225951,  0.231948,  0.237935, 
    0.243914,  0.249883,  0.255843,  0.261793,  0.267733,  0.273663,  0.279583,  0.285492, 
    0.291390,  0.297277,  0.303153,  0.309017,  0.314870,  0.320710,  0.326539,  0.332355, 
    0.338158,  0.343949,  0.349727,  0.355491,  0.361242,  0.366979,  0.372702,  0.378411, 
    0.384106,  0.389786,  0.395451,  0.401102,  0.406737,  0.412356,  0.417960,  0.423549, 
    0.429121,  0.434676,  0.440216,  0.445738,  0.451244,  0.456733,  0.462204,  0.467658, 
    0.473094,  0.478512,  0.483911,  0.489293,  0.494656,  0.500000,  0.505325,  0.510631, 
    0.515918,  0.521185,  0.526432,  0.531659,  0.536867,  0.542053,  0.547220,  0.552365, 
    0.557489,  0.562593,  0.567675,  0.572735,  0.577774,  0.582791,  0.587785,  0.592758, 
    0.597707,  0.602635,  0.607539,  0.612420,  0.617278,  0.622113,  0.626924,  0.631711, 
    0.636474,  0.641213,  0.645928,  0.650618,  0.655284,  0.659925,  0.664540,  0.669131, 
    0.673696,  0.678235,  0.682749,  0.687237,  0.691698,  0.696134,  0.700543,  0.704926, 
    0.709281,  0.713610,  0.717912,  0.722186,  0.726434,  0.730653,  0.734845,  0.739009, 
    0.743145,  0.747253,  0.751332,  0.755383,  0.759405,  0.763398,  0.767363,  0.771298, 
    0.775204,  0.779081,  0.782928,  0.786745,  0.790532,  0.794290,  0.798017,  0.801714, 
    0.805381,  0.809017,  0.812622,  0.816197,  0.819740,  0.823253,  0.826734,  0.830184, 
    0.833602,  0.836989,  0.840344,  0.843667,  0.846958,  0.850217,  0.853444,  0.856638, 
    0.859800,  0.862929,  0.866025,  0.869089,  0.872120,  0.875117,  0.878081,  0.881012, 
    0.883910,  0.886774,  0.889604,  0.892401,  0.895163,  0.897892,  0.900587,  0.903247, 
    0.905873,  0.908465,  0.911023,  0.913545,  0.916034,  0.918487,  0.920906,  0.923289, 
    0.925638,  0.927951,  0.930229,  0.932472,  0.934680,  0.936852,  0.938988,  0.941089, 
    0.943154,  0.945184,  0.947177,  0.949135,  0.951057,  0.952942,  0.954791,  0.956604, 
    0.958381,  0.960122,  0.961826,  0.963493,  0.965124,  0.966718,  0.968276,  0.969797, 
    0.971281,  0.972728,  0.974139,  0.975512,  0.976848,  0.978148,  0.979410,  0.980635, 
    0.981823,  0.982973,  0.984086,  0.985162,  0.986201,  0.987202,  0.988165,  0.989092, 
    0.989980,  0.990831,  0.991645,  0.992421,  0.993159,  0.993859,  0.994522,  0.995147, 
    0.995734,  0.996284,  0.996795,  0.997269,  0.997705,  0.998103,  0.998464,  0.998786, 
    0.999070,  0.999317,  0.999526,  0.999696,  0.999829,  0.999924,  0.999981,  1.000000, 1.000000
};

double fmodn(double x, double y)
{
    return x - y * floor(x/y);
}

double mix(double x, double y, double a) // (lerp)
{
    return x * a + y * (1.0 - a);
}

double sine(double arg)
{
    arg = fmodn(arg, 2.0 * PI); // range 0.0 to 2.0 * PI

    const int quadrant = (int)(arg / (0.5 * PI)); // range 0 to 3
    const bool flip = (quadrant > 2); // -retval
    const bool invert = (quadrant % 2); // 1.0 - arg

    arg = fmod(arg, 0.5 * PI) / (0.5 * PI); // range 0.0 to 1.0
    arg = invert ? 1.0 - arg : arg;

    const int index = arg / STEP_SIZE;
    const double rem = arg - index * STEP_SIZE;
    double retval = mix(sinLUT[index], sinLUT[index + 1], rem / STEP_SIZE);
    
    retval = flip ? -retval : retval;

    return retval;
}

#include <stdio.h>

int main()
{
    for (double i = -2.5 * PI; i < 2.5 * PI; i += PI * 0.05)
        printf("%f %fn", sin(i), sine(i));
}

Эта версия отображает arg в диапазоне от 0,0 до 1,0 перед вычислением индекса в таблице (и изменения STEP_SIZE чтобы соответствовать). Это немного медленнее, и в этом нет необходимости.

Стандартная библиотека может не гарантировать «фиксированного» времени выполнения, но это не значит, что конкретная реализация платформы не будет соответствовать вашим потребностям.

Я не слишком знаком с этим материалом; Я уверен, что есть гораздо более быстрые и точные алгоритмы …

Одно замечание заключается в том, что вы должны избавиться от всех ненужных ветвлений и повторения кода. Это плохо как для производительности, так и для обслуживания кода.

С учетом угла вы должны уметь:

• Возьмите это абсолютное значение.
• Разделить на PI / 2.
• Преобразование в беззнаковое целое с усечением десятичных знаков.

Тогда у вас либо будет индекс от 0 до 3, либо вы начнете с угла, превышающего 2 * PI (и, возможно, сначала придется разобраться с этим, если это допустимый вариант использования?).

Псевдокод:

double angle = ...;
angle /= PI/2.0;
unsigned long q = (unsigned long)angle;
// optionally some error handling/sanity checking here

quadrant(q, angle);

Если функция квадранта может использовать переменную q в качестве индекса к различным другим таблицам, содержащим множители и т. д., необходимые для ваших расчетов. Он должен сводиться к 5-6 строкам, которые у вас есть внутри каждого if-else if.

(C ++) std :: sin примерно в 9 раз быстрее, чем версия синуса LUT с плавающей запятой, закодированная @ user673679 выше. Не знаю, насколько быстрее / медленнее sinf в C.

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

Можете ли вы воспользоваться преимуществами SSE / AVX?