ここでは，DCT-I (スケーリングは除外) :
N
Ａ ＝ Σ ａ cos(πjk/N) ， 0 <= k <= N
k j=0 j
を長さ N の実数 FFT を用いて計算する方法を示します [参考文献]． その方法は多少技巧的で，まず配列 ｒ_j を
1
ｒ ＝ --- （ａ ＋ ａ ） − （ａ − ａ ） sin(πj/N)
j 2 j N-j j N-j
のようにすると，次の関係が成り立つことを利用したものです．
N-1
Ａ ＝ Σ ｒ cos(2πjk/N) ， 0 <= k <= N/2
2k j=0 j
N-1
Ａ − Ａ ＝ Σ ｒ sin(2πjk/N) ， 1 <= k <= N/2-1
2k+1 2k-1 j=0 j
したがって，配列 ａ から 配列 ｒ を計算し，長さ N の実数 FFT を 呼び出せばよいことになります．この実数 FFT の出力は，実部と虚部が Ａ の偶数項と二つの奇数項の差分になります．したがって，Ａ の奇数項に 関しては，Ａ_1 だけを別に計算しておいて，上の漸化式により逐次的に 求めなければなりません．
同様に，DST-I :
N-1
Ａ ＝ Σ ａ sin(πjk/N) ， 0 <= k <= N
k j=1 j
に関しても，配列 ｓ_j :
1
ｓ ＝ --- （ａ − ａ ） ＋ （ａ ＋ ａ ） sin(πj/N)
j 2 j N-j j N-j
を計算し，長さ N の実数 FFT を呼び出します．そして，関係式 :
N-1
Ａ ＝ Σ ｓ sin(2πjk/N) ， 1 <= k <= N/2-1
2k j=0 j
N-1
Ａ − Ａ ＝ Σ ｓ cos(2πjk/N) ， 0 <= k <= N/2
2k+1 2k-1 j=0 j
により，Ａ_k を復元します．このときの漸化式の初項は， Ａ_1 = -Ａ_-1 という関係から簡単に求められ，上の式で k=0 の項を 半分にするだけで計算できます．
この方法には，実数の FFT が直接利用できるという利点があります．逆に， 欠点として，漸化式の計算で丸め誤差が蓄積されてしまい，精度が悪くなる ということが起こります．この丸め誤差による誤差は，最悪で長さ N に 比例して大きくなります．したがって，大きい長さの計算には不向きです．
ここでは，DCT-I :
N
Ａ ＝ Σ ａ cos(πjk/N) ， 0 <= k <= N
k j=0 j
を DCT-II または DCT-III に分解する方法を示します [参考文献]． ここでは分解のみを示し，DCT-II，DCT-III の計算方法は後に示します． この分解は，Cooley-Tukey 型 FFT の分解と同じ方法であり，まず長さ N の DCT-I を N/2 の長さの DCT-I と N/2 の長さの DCT-III に分解します．
N/2
Ａ ＝ Σ （ａ ＋ ａ ） cos(πjk/(N/2))
2k j=0 j N-j
N/2-1
Ａ ＝ Σ （ａ − ａ ） cos(πj(k+1/2)/(N/2))
2k+1 j=0 j N-j
この分解を繰り返し用いて，DCT-III のみに分解します．この分解は， 周波数間引きの FFT に相当するものですが，時間間引きに相当する分解でも 同様に長さ N の DCT-I は長さ N/2 の DCT-I と長さ N/2 の DCT-II に 分解されます．このとき，DCT-I は DCT-II に分解できます．
DST-I に関しても全く同様のことが成り立ち，DST-I は DST-II または DST-III に分解できます．
まず最初に，DCT-II は，データの並べ替えだけで時間方向に 1/4 シフトした拡張 DFT にそのまま置き換えられることを示します [参考文献]．まず，DCT-II :
N-1
Ａ ＝ Σ ａ cos(π(j+1/2)k/N) ， 0 <= k < N
k j=0 j
を計算するには，入力データを次のように並べ替えます．
ｆ ＝ ａ ， 0 <= j < N/2
j 2j
ｆ ＝ ａ ， N/2 <= j < N
j 2N-2j-1
次に，時間方向に 1/4 シフトした(1/2 シフトではない)長さ N の 拡張 DFT を行います．
N-1
Ｆ ＝ Σ ｆ exp(2πi(j+1/4)k/N) ， 0 <= k <= N/2
k j=0 j
結果は
1
Ａ ＝ --- （Ｆ − i Ｆ ）
k 2 k N-k
1
Ａ ＝ ---- （Ｆ ＋ i Ｆ ）
N-k 2i k N-k
から得られます．入力データが実数のときは，-iＦ_N-k は Ｆ_k の 複素共役になるので，Ｆ_k の実部と虚部は Ａ_k と Ａ_N-k に対応します． したがって，長さ N の DCT-II の計算は，長さ N の実数の拡張 DFT の 計算そのものに置き換えられたことになります．
この方法は言われてみれば納得するのですが，何も知らないで 導出するのは困難だと思います．
次に，通常の実数 FFT を用いて拡張 DFT を計算するルーチンを示します．
| リスト2.3.2-1. 実数 FFT を用いた FCT-II |
|---|
/*
離散コサイン変換 (Unnormalized DCT-II) :
F[k]=Σ_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n),0<=k<n を計算する.
*/
#include <math.h>
void fct_ii(int n, double a[], double tmp[])
{
void rfft(int n, double a[]);
int nh, j, k;
double theta, wr, wi, xr, xi;
nh = n / 2;
for (j = 0; j < nh; j++) {
tmp[j] = a[2 * j];
tmp[n - 1 - j] = a[2 * j + 1];
}
/* ---- 1/4 Shifted DFT ---- */
rfft(n, tmp);
theta = 2 * atan(1.0) / n; /* 2*pi/(4*n) */
for (k = 1; k < nh; k++) {
wr = cos(theta * k);
wi = sin(theta * k);
xr = tmp[k];
xi = tmp[n - k];
a[k] = wr * xr - wi * xi;
a[n - k] = wr * xi + wi * xr;
}
a[0] = tmp[0];
a[nh] = tmp[nh] * cos(theta * nh);
}
|
このルーチンは，実数 FFT を呼び出し，結果に Ｗ^(j/4) を 乗じて拡張 DFT とします．また，以前に示したように 1/4 シフトした 拡張 DFT は周波数間引きの FFT の三角関数データを 1/4 ずらすことでも 計算することができます．この拡張 FFT はさらに，実数データの冗長計算を 省くことができ，DCT-II が直接 FFT で計算できることになります． 実はこの拡張 FFT を用いる方法は，対称性を利用して拡張 FFT の 計算量を 1/4 に減らす方法 (後に示す DCT/DST 専用の FFT ルーチンを用いる第一の方法) と全く同じになります．このルーチンは後に示します．
また，DCT-III は DCT-II の逆変換なので，単に計算手順を 逆にすればよく，DCT-III :
N-1
Ａ ＝ Σ ａ cos(πj(k+1/2)/N) ， 0 <= k < N
k j=0 j
は，実数の拡張 DFT の逆変換 :
N/2 N/2-1
Ｆ ＝ Σ ａ cos(2πj(k+1/4)/N) ＋ Σ ａ sin(2πj(k+1/4)/N)
k j=0 j j=1 N-j
になります．結果は
Ａ ＝ Ｆ ， 0 <= k < N/2
2k k
Ａ ＝ Ｆ ， N/2 <= k < N
2N-2k-1 k
となります．実数の FFT の逆変換を用いたルーチンは次のようになります．
| リスト2.3.2-2. 実数 FFT の逆変換を用いた FCT-III |
|---|
/*
逆離散コサイン変換 (Unnormalized DCT-III) :
F[k]=Σ_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n),0<=k<n を計算する.
*/
#include <math.h>
void fct_iii(int n, double a[], double tmp[])
{
void irfft(int n, double a[]);
int nh, j, k;
double theta, wr, wi, xr, xi;
nh = n / 2;
theta = -2 * atan(1.0) / n; /* 2*pi/(4*n) */
for (k = 1; k < nh; k++) {
wr = cos(theta * k);
wi = sin(theta * k);
xr = a[k];
xi = a[n - k];
tmp[k] = wr * xr - wi * xi;
tmp[n - k] = wr * xi + wi * xr;
}
tmp[0] = a[0] * 2;
tmp[nh] = a[nh] * 2 * cos(theta * nh);
irfft(n, tmp);
for (j = 0; j < nh; j++) {
a[2 * j] = tmp[j];
a[2 * j + 1] = tmp[n - 1 - j];
}
}
|
さらにこの方法は，多次元に拡張することができ，多次元 DCT-II, DCT-III は多次元の 1/4 シフトした拡張 DFT の計算に置き換えられます [参考文献]．
ここでは，DCT-II, DCT-III を長さ N の実数 FFT を用いて計算する 別の方法を示します．この方法は多少曲芸的ですが，導出方法は第一の 方法より簡単です．導出の概略は，DCT-III の 1/2 だけずれた項を展開して 半分の長さの DCT-I と DST-I のペアに分解し，さらに対称性を使って 実数データの DFT に置き換えるというものです．結果だけ示すと，DCT-III :
N-1
Ａ ＝ Σ ａ cos(πj(k+1/2)/N) ， 0 <= k < N
k j=0 j
は，入力データに対する演算 :
1+i j
ｆ ＝ Re（(ａ − i ａ ) ----- exp(2πi----)）， 1 <= j <= N/2
j j N-j 2 4N
1+i j
ｆ ＝ Im（(ａ − i ａ ) ----- exp(2πi----)）， 1 <= j < N/2
N-j j N-j 2 4N
ｆ ＝ ａ
0 0
を行うことで，通常の実数データの FFT に変換されます．配列 ｆ の 実数 DFT :
N-1
Ｆ ＝ Σ ｆ exp(2πijk/N) ， 0 <= k <= N/2
k j=0 j
を計算すると，DCT-III の結果は
Ａ ＝ Re（ (1-i) Ｆ ） ， 0 <= k < N/2
2k k
Ａ ＝ -Im（ (1-i) Ｆ ） ， 0 < k <= N/2
2k-1 k
で復元されます．ちなみに，この式の確認は難しくはないのですが， 式を全部展開してみないとわからないと思います．
この方法の演算量は，第一の方法に比べて実数加算が N-2 回だけ 多くなります．しかし，実数データの FFT の出力の実部と虚部が隣り合わせに なるならば，この方法はデータの並べ替えが不用になり，in-place 演算になる という利点があります．一般に，データの並べ替えは単純な計算よりも 時間のかかる場合が多いので，この方法は多くの状況で有利になるようです．
前と同様に，DCT-II に対する方法は，この手続きを逆にすることで 得られます．
ここでは，DCT-IV を半分の長さの DCT-II または DCT-III に分解して 計算する方法を示します[参考文献]． この方法は，DCT-IV :
N-1
Ａ ＝ Σ ａ cos(π(j+1/2)(k+1/2)/N) ， 0 <= k < N
k j=0 j
が次のように変形できることを用います．
1 N-1
---（Ａ ＋ Ａ ）＝ Σ ａ cos(π(j+1/2)/2N) cos(π(j+1/2)2k/N)
2 2k-1 2k j=0 j
1 N-1
---（Ａ − Ａ ）＝ Σ ａ sin(π(j+1/2)/2N) sin(π(j+1/2)2k/N)
2 2k-1 2k j=0 j
したがって，DCT-IV は，半分の長さの DCT-II と DST-II :
1 N/2-1
---（Ａ ＋ Ａ ）＝ Σ ｒ cos(π(j+1/2)k/(N/2))
2 2k-1 2k j=0 j
1 N/2-1
---（Ａ − Ａ ）＝ Σ ｓ sin(π(j+1/2)k/(N/2))
2 2k-1 2k j=0 j
に分解することができます．ここで，配列 ｒ と ｓ は
ｒ ＝ ａ cos(π(j+1/2)/2N) ＋ ａ sin(π(j+1/2)/2N)
j j N-1-j
ｓ ＝ ａ sin(π(j+1/2)/2N) − ａ cos(π(j+1/2)/2N)
j j N-1-j
です．この演算は複素数乗算に相当する演算になっていることに 注意してください．さらに，DST-II は入力データの符号の付け替えと 出力データの並べ替えだけで DCT-II :
1 N/2-1 j
---（Ａ − Ａ ）＝ Σ ｓ (-1) cos(π(j+1/2)(N/2-k)/(N/2))
2 2k-1 2k j=0 j
になります．
前と同様に，DCT-III に分解する方法は，この手続きを逆にすることで 得られます．