Z変換から導いた漸化式で正弦波を高速に生成

マイコンで三角関数を計算するのは重い処理ですが、次の漸化式を用いると離散正弦波の数列を高速に生成することができます。漸化式なので生成を繰り返すと誤差が蓄積することに注意。

離散正弦波の数列 x(n) は
正弦波の周波数を f 、サンプリング周波数を k とすると
x(0) = 0
x(1) = sin( 2πf / k )
n≧2 のとき
x(n) = a x(n-1) - x(n-2)
ただし、a = 2 cos( 2πf / k )

下記のページに詳しい説明があります。

• 漸化式でsin波を高速に生成する方法 - nabeの雑記帳

違う見方をすると、正弦波のZ変換からこの漸化式を導くこともできます。
(たぶん本質的には同じことです。)

• メモ１　メモ２

実際に計算した結果のグラフです。

計算に用いたプログラムです。

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

const float PI = 3.14159265359;

// Z変換から導いた漸化式で離散正弦波の数列を生成するクラス
class ZtrSin
{
public:
    float xz1;  // x[n-1] : 前回値
    float xz2;  // x[n-2] : 前々回値
    float a;    // 漸化式の係数

public:
    // 初期化 (漸化式の初期値と係数の計算)
    // f : 正弦波の周波数[Hz]
    // k : サンプリング周波数[Hz]
    // phi : 初期位相[rad]
    ZtrSin(float f, float k, float phi)
    {
        xz2 = sinf(phi);
        xz1 = sinf(phi + 2 * PI * f / k);
        a = 2 * cos(2 * PI * f / k);
    }

    // 次の値を生成 (漸化式の計算)
    float get() {
        float x = a * xz1 - xz2;
        xz2 = xz1;
        xz1 = x;
        return x;
    }
};

int main(void)
{
    // f=10Hz, k=320Hz, phi=0の離散正弦波の数列を2周期ぶん出力 
    ZtrSin ztrsin(10, 320, 0);
    printf("%f\n", ztrsin.xz2);
    printf("%f\n", ztrsin.xz1);
    for (int n = 2; n <= 64; n++) {
        printf("%f\n", ztrsin.get());
    }

    return 0;
}