FFT向け三角関数表の作成

以下は奥村晴彦著 "C言語による最新アルゴリズム事典"に記載のあるFFT計算のコードの一部です。
FFTを実行する際に利用する三角関数表(sinテーブル)を作成しています。
sinを一度求めたら、それを使って次々にcos,sinのn倍角の値を計算しています。
が、書籍中に特に説明はなく、なぜこれで計算できるのかすぐにはわかりません。
（FFTのソースコードはここにあります。）

/*
  関数{\tt fft()}の下請けとして三角関数表を作る.
*/
static void make_sintbl(int n, double sintbl[])
{
    int i, n2, n4, n8;
    double c, s, dc, ds, t;

    n2 = n / 2;  n4 = n / 4;  n8 = n / 8;
    t = sin(PI / n);
    dc = 2 * t * t;  ds = sqrt(dc * (2 - dc));
    t = 2 * dc;  c = sintbl[n4] = 1;  s = sintbl[0] = 0;
    for (i = 1; i < n8; i++) {
        c -= dc;  dc += t * c;
        s += ds;  ds -= t * s;
        sintbl[i] = s;  sintbl[n4 - i] = c;
    }
    if (n8 != 0) sintbl[n8] = sqrt(0.5);
    for (i = 0; i < n4; i++) {
        sintbl[n2 - i] = sintbl[i];
        sintbl[i + n2] = - sintbl[i];
    }
}

三角関数を計算している中心的な部分

上のコードの最も肝心なところを以下に抜粋しました。

    t = sin(PI / n);
    dc = 2 * t * t;  ds = sqrt(dc * (2 - dc));
    t = 2 * dc;  c = sintbl[n4] = 1;  s = sintbl[0] = 0;
    for (i = 1; i < n8; i++) {
        c -= dc;  dc += t * c; // cosに関する計算
        s += ds;  ds -= t * s; // sinに関する計算
        sintbl[i] = s;  sintbl[n4 - i] = c; // 求めたcos,sinの値をテーブルへ設定
    }

上の部分で$\theta =2\pi/N $ として
変数 c に$\cos(n\theta)$、変数 s に$\sin(n\theta)$ が計算されているようです。（実際計算できている）
これについて以下で示しました。

チェビシェフ多項式

$\cos(\theta)$ のn倍角の値はチェビシェフ多項式で表せる。

\begin{align}
\cos(1t) &= \cos(t) \\
\cos(2t) &= \cos(t)\cos(t) - \sin(t)\sin(t)=2(\cos(t))^2 - 1\\
\cos(3t) &= 4(\cos(t))^3-3\cos(t) \\
\vdots&
\end{align}

以下の式で$x = \cos(t)$ とすれば、

\begin{align}
T_0(x)& = 1 \\
T_1(x)& = x \\
T_2(x)& = 2x^2-1 \\
T_3(x)& = 4x^3-3x \\
\vdots&
\end{align}

$T_n(x)$を漸化式で表せば

T_{n+1}(x) = 2xT_n(x)-T_{n-1}(x)

と表せる。これを利用する。

cosのn倍角について

$k = 2\sin^2(\theta/2)$とすると
$\cos(\theta) = 1 - 2\sin^2(\theta/2) = 1 - k$
また、
$ T_n(\cos\theta) $ を $T_n$と書けばチェビシェフ多項式を使って

\begin{align}
T_0 &= 1 \\
T_1 &= 1 - k = \cos(\theta) \\
T_2 &= 2(1 - k)T_1 - T_0 \\
    &= 2T_1 - 2kT_1 - T_0 = T_1 - 2kT_1 + T_1 - T_0 \\
    &= T_1 - 2kT_1 - k \\
T_3 &= 2(1 - k)T_2 - T_1 \\
    &= 2T_2 - 2kT_2 - T_1 = T_2 - 2kT_2 + T_2 - T_1 \\
    &= T_2 - 2kT_2 - 2kT_1 - k \\
\vdots&
\end{align}

以上のように書ける。これを見ると以下の式が成り立っているように見える。

T_{n+1} = T_n - \sum_{i=1}^{n} 2kT_i - k  \tag{1}

式(1)がn>=1で成り立つことの証明

式$(1)$は、$n=1$のとき成り立つ。
$n=p$のとき成り立つとすると、$n=p+1$は、

\begin{align}
T_{p+2} &= 2(1 - k)T_{p+1} - T_p \\
        &= 2T_{p+1} - 2kT_{p+1} - T_p \\
        &= T_{p+1} - 2kT_{p+1} + T_{p+1} - T_{p} \\
        &= T_{p+1} - 2kT_{p+1} + (T_p - \sum_{i=1}^{p} 2kT_i - k) - T_p \\
        &= T_{p+1} - \sum_{i=1}^{p+1} 2kT_i - k
\end{align}

つまり、$n=p+1$で成り立つ。

プログラム中のcos計算との関係

c -= dc;  dc += t * c;

T_{n+1} = T_n - \sum_{i=1}^{n} 2kT_i - k  \tag{1}

上記のコードの、変数 c は式$(1)$に従って計算されていて、$T_n$に相当する。
変数 t は$2k$に相当し、変数dcは

\sum_{i=1}^{n} 2kT_i + k

に相当することがわかる。

sinのn倍角について

sinについてもcosと同様に求める。
$k = 2\sin^2(\theta/2) $とし、
$l = \sqrt{k(2 - k)}$
とすると、

\begin{align}
l &= \sqrt{k(2 - k)} = \sqrt{2\sin^2(\theta/2)(2-2\sin^2(\theta/2)} \\
&=\sqrt{4\sin^2(\theta/2)-4\sin^4(\theta/2)} \\
&=2\sin(\theta/2)\cos(\theta/2) \\
&=\sin(\theta) \\
&=\cos(\pi/2-\theta)
\end{align}

ここで、チェビシェフ多項式$T_n$を使って$S_n$を

S_n = T_{N/4-n}

として、$\pi/2$から右回りに$\theta$角ごとの$\cos$を求める。
これはつまり、$\sin(n\theta)$を求めることになる。

\begin{align}
S_0 &= \cos(\pi/2) = 0 \\
S_1 &= \cos(\pi/2 - \theta) = \sin(\theta) = l \\
S_2 &= 2(1 - k)S_1 - S_0 \\
&= 2S_1 - 2kS_1 - S_0 = S_1-2kS_1+S_1-S_0 \\
&= S_1 - 2kS_1 + l \\
S_3&= 2(1-k)S_2-S1\\
&=2S_2-2kS_2-S_1 =S_2-2kS_2+S_2-S_1  \\
&= S_2-2kS_2-2kS_1+l \\
\vdots&
\end{align}

これを見ると以下の式が成り立っているように見える。

S_{n+1} = S_n + \sum_{i=1}^{n} -2kS_i + l  \tag{2}

式(2)がn>=1で成り立つことの証明

式$(2)$は$n=1$のときに成り立つ。$n=p$のときに成り立つとすると、
$n=p+1$について

\begin{align}
S_{p+2} &= 2(1 - k)S_{p+1} - S_p \\
        &= 2S_{p+1} - 2kS_{p+1} - S_p \\
        &= S_{p+1} - 2kS_{p+1} + S_{p+1} - S_{p} \\
        &= S_{p+1} - 2kS_{p+1} + (S_p + \sum_{i=1}^{p} -2kS_i + l) - S_p \\
        &= S_{p+1} + \sum_{i=1}^{p+1} -2kS_i + l
\end{align}

つまり、$n=p+1$で成り立つ。

プログラム中のsin計算との関係

s += ds;  ds -= t * s;

S_{n+1} = S_n + \sum_{i=1}^{n} -2kS_i + l  \tag{2}

上記のコードの、変数 s は式$(2)$に従って計算されていて、$S_n$に相当する。
変数 t は$2k$に相当し、変数 ds は式$(2)$で考えれば、

\sum_{i=1}^{n} -2kS_i + l

に相当することがわかる。

以上で三角関数テーブルの計算方法がわかった。