Как вычислить синусоиду с точностью по времени

Question

случая использования для генерации синусоиды для цифрового синтеза, поэтому нам необходимо вычислить все значения sin (дт), где:

т является целым числом номер, представляющий номер выборки. Это переменная. Диапазон составляет от 0 до 158 760 000 за один час качества CD.

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

Цель заключается в обеспечении высокой точности с традиционными INT и двойных типов данных. Производительность не важна.

Наивная реализация:

double next() 
{ 
    t++; 
    return sin(((double) t) * (d)); 
} 

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

Усовершенствованный вариант заключается в следующем:

double next() 
{ 
    d_sum += d; 
    if (d_sum >= (M_PI*2)) d_sum -= (M_PI*2); 

    return sin(d_sum); 
} 

Здесь я убеждаюсь предоставлять номера в диапазоне от 0 до 2 * пи функции «грех».

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

Вопрос в том, как повысить точность.

---

Приложение 1

"точность получает снижается из-за больших чисел, предоставленных "" функции" грех:

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

#define TEST  (300000006.7846112) 
#define TEST_MOD (0.0463259891528704262050786960234519968548937998410258872449766) 
#define SIN_TEST (0.0463094209176730795999323058165987662490610492247070175523420) 

int main() 
{ 
    double a = sin(TEST); 
    double b = sin(TEST_MOD); 

    printf("a=%0.20f \n" , a); 
    printf("diff=%0.20f \n" , a - SIN_TEST); 
    printf("b=%0.20f \n" , b); 
    printf("diff=%0.20f \n" , b - SIN_TEST); 
    return 0; 
} 

Выход:

a=0.04630944601888796475 
diff=0.00000002510121488442 
b=0.04630942091767308033 
diff=0.00000000000000000000 

Answer 1

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

Sin(A + d) = Sin(A) * Cos(d) + Cos(A) * Sin(d) 

Здесь мы должны хранить и значение косинуса обновления слишком и хранить константы (при заданном дельта) факторы Cos (d) и Sin (г).

Теперь о точности: косинус (d) для малых d очень близок к 1, поэтому существует риск потери точности (в цифрах всего лишь несколько значащих цифр, например 0.99999987).Чтобы преодолеть эту проблему, мы можем хранить постоянные факторы, как

dc = Cos(d) - 1 = - 2 * Sin(d/2)^2 
ds = Sin(d) 

использования других формул, чтобы обновить текущее значение
(здесь sa = Sin(A) для текущего значения, ca = Cos(A) для текущего значения)

ts = sa //remember last values 
tc = ca 
sa = sa * dc + ca * ds 
ca = ca * dc - ts * ds 
sa = sa + ts 
ca = ca + tc 

P.S. Некоторые реализации БПФ периодически (каждый шаг К) обновляют значения sa и ca через триггер. чтобы избежать накопления ошибок.

Пример результата. Расчеты в парном разряде.

d=0.000125 
800000000 iterations 
finish angle 100000 radians 

          cos    sin 
described method  -0.99936080743598 0.03574879796994 
Cos,Sin(100000)   -0.99936080743821 0.03574879797202 
windows Calc   -0.9993608074382124518911354141448 
          0.03574879797201650931647050069581   

Answer 2

Sin (х) = Sin (х + 2N ∙ π), так что проблема может быть сводилась точно находить небольшое число, которое равно большое число х по модулю 2π.

Например, -1,61059759 ≅ 256 мод 2π, и вы можете рассчитывать sin(-1.61059759) с большей точностью, чем sin(256)

Так давайте выберем некоторое целое число, чтобы работать с, 256. Сначала найти небольшие числа, которые равны полномочия 256, по модулю 2π:

// to be calculated once for a given frequency 
// approximate hard-coded numbers for d = 1 below: 
double modB = -1.61059759; // = 256 mod (2π/d) 
double modC = 2.37724612; // = 256² mod (2π/d) 
double modD = -0.89396887; // = 256³ mod (2π/d) 

, а затем разделить свой индекс как число в базе 256:

// split into a base 256 representation 
int a = i   & 0xff; 
int b = (i >> 8) & 0xff; 
int c = (i >> 16) & 0xff; 
int d = (i >> 24) & 0xff; 

Теперь вы можете найти гораздо меньшее число x, равное i по модулю 2π/д

// use our smaller constants instead of the powers of 256 
double x = a + modB * b + modC * c + modD * d; 
double the_answer = sin(d * x); 

Для различных значений д вы должны рассчитать различные значения modB, modC и modD, которые равны тем мощностям 256, но по модулю (2π/d). Вы можете использовать библиотеку высокой точности для этих двух вычислений.

Answer 3

Шкала вверх период до 2^64, и делать умножение с использованием целочисленной арифметики:

// constants: 
double uint64Max = pow(2.0, 64.0); 
double sinFactor = 2 * M_PI/(uint64Max); 

// scale the period of the waveform up to 2^64 
uint64_t multiplier = (uint64_t) floor(0.5 + uint64Max * d/(2.0 * M_PI)); 

// multiplication with index (implicitly modulo 2^64) 
uint64_t x = i * multiplier; 

// scale 2^64 down to 2π 
double value = sin((double)x * sinFactor); 

Пока ваш период не миллиарды образцов, точность multiplier будет достаточно хорошо.

Answer 4

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

double next() { 
    t0 += 1.0; 
    d_sum = t0 * d; 
    if (d_sum > 2.0 * M_PI) { 
     t0 -= ((2.0 * M_PI)/d); 
    } 
    return (sin(d_sum)); 
} 

Answer 5

Для гипер точности, ОП имеет 2 проблемы:

1. умножение d по n и поддержанию большой точности, чем double.На это отвечает первая часть ниже.

2. Выполнение операции mod периода. Простое решение состоит в том, чтобы использовать градусы, а затем mod 360, достаточно легко сделать точно. Для 2*π больших углов сложно, поскольку это требует значение 2*π с около 27 более битов точности, чем (double) 2.0 * M_PI

---

Используйте 2 double с представлять d.

Предположим, что 32-разрядные int и binary64double. Таким образом, double имеет точность 53 бит.

0 <= n <= 158,760,000 который составляет около 2 ^27.2. Так как double может обрабатывать 53-битные целые без знака непрерывно и точно, 53-28 -> 25, любой double с 25 значащими битами может быть умножен на n и по-прежнему быть точным.

Сегмент d на 2 double с dmsb,dlsb, 25-значные цифры и 28 наименее.

int exp; 
double dmsb = frexp(d, &exp); // exact result 
dmsb = floor(dmsb * POW2_25); // exact result 
dmsb /= POW2_25;    // exact result 
dmsb *= pow(2, exp);   // exact result 
double dlsb = d - dmsb;  // exact result 

Тогда каждое умножение (или последовательное добавление) из dmsb*n будет точным. (это важная часть.) dlsb*n будет ошибкой только в наименьших количествах бит.

double next() 
{ 
    d_sum_msb += dmsb; // exact 
    d_sum_lsb += dlsb; 
    double angle = fmod(d_sum_msb, M_PI*2); // exact 
    angle += fmod(d_sum_lsb, M_PI*2); 
    return sin(angle); 
} 

Примечание: fmod(x,y) результаты, как ожидается, будет точной дать точную x,y.

---

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

#define AS_n 158760000 
double AS_d = 300000006.7846112/AS_n; 
double AS_d_sum_msb = 0.0; 
double AS_d_sum_lsb = 0.0; 
double AS_dmsb = 0.0; 
double AS_dlsb = 0.0; 

double next() { 
    AS_d_sum_msb += AS_dmsb; // exact 
    AS_d_sum_lsb += AS_dlsb; 
    double angle = fmod(AS_d_sum_msb, M_PI * 2); // exact 
    angle += fmod(AS_d_sum_lsb, M_PI * 2); 
    return sin(angle); 
} 

#define POW2_25 (1U << 25) 

int main(void) { 
    int exp; 
    AS_dmsb = frexp(AS_d, &exp);   // exact result 
    AS_dmsb = floor(AS_dmsb * POW2_25); // exact result 
    AS_dmsb /= POW2_25;     // exact result 
    AS_dmsb *= pow(2, exp);    // exact result 
    AS_dlsb = AS_d - AS_dmsb;   // exact result 

    double y; 
    for (long i = 0; i < AS_n; i++) 
    y = next(); 
    printf("%.20f\n", y); 
} 

Выходные

0.04630942695385031893 

---

Использование градусов

Рекомендуется использовать в качестве градусов 360 градусов по точные периода и M_PI*2 радианах приложение roximation. C не может представлять π точно.

Если ОП все еще хочет использовать радианы, для дальнейшего понимания о выполнении мод в я, см Good to the Last Bit