7

使用例は、デジタル合成用の正弦波を生成することです。そのため、sin(dt) のすべての値を計算する必要があります。

tはサンプル番号を表す整数です。これは可変です。CD 品質の 1 時間のサウンドの範囲は 0 ～ 158,760,000 です。

dは double で、角度のデルタを表します。これは一定です。範囲は次のとおりです: 0 より大きい、pi より小さい。

目標は、従来の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);
}

ここでは、「sin」関数に 0 から 2*pi の範囲の数値を指定するようにします。

しかし、問題は、dが小さい場合、小さな追加が多くなり、毎回精度が低下することです。

ここで問題となるのは、精度を向上させる方法です。

付録1

「「sin」関数に大きな数値が提供されるため、精度が低下します」:

#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
4

5 に答える 5

1

周期を 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十分です。

于 2016-06-08T10:05:54.230 に答える
0

超精度の場合、OP には 2 つの問題があります。

  1. を掛けdてn、より高い精度を維持しますdouble。それは以下の最初の部分で答えられます。

  2. 期間のを実行しmodます。簡単な解決策は、度を使用してから を使用mod 360することです。正確に行うのは簡単です。大きな角度を行う2*πには注意が必要2*πです。(double) 2.0 * M_PI

doubleを表すには2 を使用しdます。

int32 ビットとbinary64 を想定してみましょうdouble。53doubleビットの精度もあります。

0 <= n <= 158,760,000これは約 2 27.2です。double53 ビットの符号なし整数を継続的かつ正確に処理できるため、53-28 --> 25 のように、有効doubleビット数が 25 しかないものを乗算しnても正確です。

最上位 25 桁と最下位 28 桁のd2 doublesにセグメント化します。dmsb,dlsb

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 がまだラジアンを使用したい場合、π の mod の実行に関する詳細な洞察については、最後のビットへのグッドを参照してください。

于 2016-06-13T18:54:47.183 に答える