Как рассчитать синусоидальную волну с точностью до времени
Вариант использования - генерация синусоидальной волны для цифрового синтеза, поэтому нам нужно вычислить все значения sin(d t), где:
t является целым числом, представляющим номер образца. Это переменная. Диапазон от 0 до 158 760000 за один час звучания качества CD.
d является двойным, представляя дельту угла. Это постоянно. И диапазон: больше 0, меньше чем пи.
Цель состоит в том, чтобы достичь высокой точности с традиционными типами данных int и double. Производительность не важна.
Наивная реализация это:
double next()
{
t++;
return sin( ((double) t) * (d) );
}
Но проблема в том, что когда t увеличивается, точность уменьшается, потому что большие числа предоставлены функции "sin".
Улучшенная версия:
double next()
{
d_sum += d;
if (d_sum >= (M_PI*2)) d_sum -= (M_PI*2);
return sin(d_sum);
}
Здесь я обязательно предоставляю числа в диапазоне от 0 до 2*pi функции "sin".
Но теперь проблема в том, что когда d мало, есть много небольших дополнений, которые каждый раз снижают точность.
Вопрос здесь в том, как повысить точность.
Приложение 1
msgstr "точность снижается, потому что большие числа предоставлены для функции" греха "
#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
5 ответов
Вы можете попробовать подход, который используется в некоторых реализациях быстрого преобразования Фурье. Значения тригонометрической функции рассчитываются на основе предыдущих значений и дельты.
Sin(A + d) = Sin(A) * Cos(d) + Cos(A) * Sin(d)
Здесь мы также должны сохранить и обновить значение косинуса и сохранить постоянные (для заданной дельты) коэффициенты Cos(d) и Sin(d).
Теперь о точности: косинус (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
PS Некоторые реализации FFT периодически (каждые K шагов) обновляются 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
sin (x) = sin (x + 2N ∙ π), поэтому задачу можно свести к точному поиску небольшого числа, равного большому числу x по модулю 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π/d
// 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);
Для разных значений d вам придется рассчитывать разные значения modB, modC а также modD, которые равны этим степеням 256, но по модулю (2π / d). Вы можете использовать высокоточную библиотеку для этих двух вычислений.
Увеличьте период до 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 будет достаточно хорошо
Для повышенной точности у OP есть 2 проблемы:
умножения
dотnи поддерживая большую точность, чемdouble, Это ответ в первой части ниже.Выполнение
modпериода. Простое решение заключается в использовании градусов, а затемmod 360, достаточно легко сделать точно. Сделать2*πбольших углов сложно, так как для этого нужно значение2*πс точностью примерно на 27 бит больше, чем(double) 2.0 * M_PI
Используйте 2 double с представлять d ,
Допустим, 32-битный int и двоичный код64 double, Так double имеет 53-битную точность.
0 <= n <= 158,760,000 что составляет около 2 27,2. поскольку double может обрабатывать 53-разрядные целые числа без знака непрерывно и точно, 53-28 -> 25, любой double только с 25 значащими битами можно умножить на n и все еще быть точным.
сегмент d в 2 double s 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 радианы это приближение. C не может точно представлять π.
Если OP по-прежнему хочет использовать радианы, для получения дополнительной информации о выполнении мода π см. Good to the Last Bit
Следующий код удерживает входные данные для функции 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));
}