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
に相当することがわかる。
以上で三角関数テーブルの計算方法がわかった。