ディジタルロックインアンプに必要な計算精度

FPGAでロックインアンプを作っているのですが、同じ振幅でも位相を変えていくと、計算される振幅が10uVくらい変わるという現象が起きています。

問題を簡単にするために、C言語で同じ処理を書いてCPUで実行してみることにしました。

const int PERIOD = 1280;
const int MAXITER = 1000;
const double A = 1;
for (double phase = 0; phase < 360; phase += 1)
{
    double sum_sin = 0;
    double sum_cos = 0; // ①
    for (int iter = 0; iter < MAXITER; iter++)
    {
        for (int p = 0; p < PERIOD; p++)
        {
            double sig = (A * sin(((double)p / PERIOD + phase / 360) * 2 * PI)) * 32768; // ②
            double S = sin((double)p / PERIOD * 2 * PI) * 32768;
            double C = cos((double)p / PERIOD * 2 * PI) * 32768; // ③
            sum_sin += ((long long)S * sig / 32768);
            sum_cos += ((long long)C * sig / 32768);
        }
    }
    double X = sum_sin / 32768. / PERIOD / MAXITER * 2;
    double Y = sum_cos / 32768. / PERIOD / MAXITER * 2;
    printf("%f\t", phase);
    printf("%f\t", (sqrt(X * X + Y * Y) - A) * 1000);
}

上のプログラムは正弦波と正弦波を掛け合わせて、一周期(PERIOD)の間積分するという計算をITER回繰り返すものです。それをphase=0~359まで繰り返します。

①②③の行にあるdoubleや、long longという型が結果にどのような影響を出すかを調べます。

まず、すべてがdouble型で計算した場合の結果です。見た感じでは、誤差はほぼゼロです。

拡大すると演算の誤差は概ね10-9レベルであることがわかります。

次に③の部分をlongにした場合。つまり、SIN、COSを固定小数点32bit型にした場合です。

-20uVの誤差が出ていますが、周期的構造は見当たりません。

次に②のsigをlong型にした場合。周期的な構造が見えてきました。

①のsum用の変数をlong longにした場合。

逆に、①のsumはlong long型だけれども、SIN、COSとsigをdoubleにした場合。

やはり、周期的な誤差が出てきます。

①と③、つまりsumとSIN,COSはdoubleで、sigはlongの場合。

様子が変わってきますが、より細かい周期のギザギザとなります。

最後は①②③のすべてをlong longにした場合ですが、32768倍ではなく2147483648倍して、32bitの固定小数点として計算してみます。すると、

ほぼ、ゼロでした。

つまり、固定小数点で演算するとしても32bit×32bit=64bitくらいの精度で計算すれば誤差は出てこないと思われます。