実数データの DCT を計算するもっとも単純な方法は，データの長さを 倍にしてデータを対称にしてから，実数データに対する拡張 FFT で計算する ことです．しかし，これは半分の無駄な計算が含まれます．ここでは， データの対称性を用いて，演算量と記憶領域を実数データの拡張 FFT の 1/2 (複素数データの拡張 FFT の1/4) にする方法について示します [参考文献]．
最初に，DCT-II :
N-1
Ａ ＝ Σ ａ cos(π(j+1/2)k/N) ， 0 <= k < N
k j=0 j
は，対称なデータの拡張 DFT :
2N-1
Ａ ＝ Σ ｂ exp(2π(j+1/2)k/2N) ， 0 <= k < N
k j=0 j
に相当します．ここで，配列 ｂ は ａ を対称にしたもので，
1
ｂ ＝ ｂ ＝ --- ａ ， 0 <= j < N
j 2N-1-j 2 j
です．この対称性は，実数データに対する拡張 FFT の log(N) 回ある すべての段のデータに複素共役対称として波及していきます．そして， 最後の段で Ａ_k が Ａ_k 自身の複素共役対称となり，実数に戻ります． この冗長な演算と記憶領域は， 実数データの拡張 FFTの リスト 2.2-2 の拡張 FFT のバタフライループの長さを半分にすることで簡単に 取り除くことができます．
次に，対称性を用いて実数に対する拡張 FFT の演算量と記憶領域を 半分にした DCT-II のルーチンを示します．
| リスト2.3.3-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[])
{
int m, mh, mq, i, j, k, irev, jr, ji, kr, ki;
double theta, wr, wi, xr, xi;
/* ---- scrambler ---- */
i = 0;
for (j = 1; j < n - 1; j++) {
for (k = n >> 1; k > (i ^= k); k >>= 1);
if (j < i) {
xr = a[j];
a[j] = a[i];
a[i] = xr;
}
}
theta = -4 * atan(1.0) / n;
for (mh = 2; (m = mh << 1) <= n; mh = m) {
mq = mh >> 1;
/* ---- real to real & complex butterflies ---- */
irev = 0;
for (jr = 0; jr < n; jr += m) {
wr = cos(0.5 * theta * (irev + mq));
wi = sin(0.5 * theta * (irev + mq));
for (k = n >> 2; k > (irev ^= k); k >>= 1);
ki = jr + mq;
kr = n - ki;
ji = kr - mq;
xr = a[jr] - a[kr];
xi = -a[ji] + a[ki];
a[jr] += a[kr];
a[ji] += a[ki];
a[ki] = wr * xr + wi * xi;
a[kr] = wr * xi - wi * xr;
}
if (mh == 2) continue;
/* ---- complex to complex butterflies ---- */
irev = 0;
for (i = 0; i < n; i += m) {
wr = cos(theta * (irev + mq));
wi = sin(theta * (irev + mq));
for (k = n >> 2; k > (irev ^= k); k >>= 1);
for (j = 1; j < mq; j++) {
jr = i + j;
ki = i + mh - j;
ji = n - jr;
kr = n - ki;
xr = a[jr] - a[kr];
xi = -a[ji] - a[ki];
a[jr] += a[kr];
a[ji] -= a[ki];
a[ki] = wr * xr + wi * xi;
a[kr] = wr * xi - wi * xr;
}
}
}
if (n > 1) {
m = n >> 1;
xr = a[0] - a[m];
a[0] += a[m];
a[m] = xr * cos(0.5 * theta * m);
}
}
|
このルーチンの N=8 のときのデータフローは次のようになります．
|
| 図2.3.3-1. 拡張 FFT の冗長な演算を取り除いた FCT-II |
|---|
次に，この手続きを逆にした DCT-III に対するルーチンを示します．
| リスト2.3.3-2. 拡張 FFT の冗長な演算を取り除いた FCT-III |
|---|
/*
逆離散コサイン変換 (Unnormalized DCT-III) :
a[k]=Σ_j=0^n-1 F[j]*cos(pi*j*(k+1/2)/n),0<=k<n を計算する.
*/
#include <math.h>
void fct_iii(int n, double a[])
{
int m, mh, mq, i, j, k, irev, jr, ji, kr, ki;
double theta, wr, wi, xr, xi;
theta = 4 * atan(1.0) / n;
if (n > 1) {
m = n >> 1;
xr = a[m] * cos(0.5 * theta * m);
a[m] = a[0] - xr;
a[0] += xr;
}
for (m = n; (mh = m >> 1) >= 2; m = mh) {
mq = mh >> 1;
/* ---- real & complex to real butterflies ---- */
irev = 0;
for (jr = 0; jr < n; jr += m) {
wr = cos(0.5 * theta * (irev + mq));
wi = sin(0.5 * theta * (irev + mq));
for (k = n >> 2; k > (irev ^= k); k >>= 1);
ki = jr + mq;
kr = n - ki;
ji = kr - mq;
xr = wr * a[ki] + wi * a[kr];
xi = wr * a[kr] - wi * a[ki];
a[kr] = a[jr] - xr;
a[ki] = a[ji] + xi;
a[jr] += xr;
a[ji] -= xi;
}
if (mh == 2) continue;
/* ---- complex to complex butterflies ---- */
irev = 0;
for (i = 0; i < n; i += m) {
wr = cos(theta * (irev + mq));
wi = sin(theta * (irev + mq));
for (k = n >> 2; k > (irev ^= k); k >>= 1);
for (j = 1; j < mq; j++) {
jr = i + j;
ki = i + mh - j;
ji = n - jr;
kr = n - ki;
xr = wr * a[ki] + wi * a[kr];
xi = wr * a[kr] - wi * a[ki];
a[kr] = a[jr] - xr;
a[ki] = -a[ji] - xi;
a[jr] += xr;
a[ji] -= xi;
}
}
}
/* ---- unscrambler ---- */
i = 0;
for (j = 1; j < n - 1; j++) {
for (k = n >> 1; k > (i ^= k); k >>= 1);
if (j < i) {
xr = a[j];
a[j] = a[i];
a[i] = xr;
}
}
}
|
この方法は，基数を大きくしたり Split-Radix FFT にしたりすることで， 演算量はさらに減らすことができます．さらに， 混合基数にすることで任意の 長さの DCT が計算できます．
ここでは，Cooley-Tukey 型 FFT と同様の考え方で，長さ N の DCT-II を 二つの N/2 の長さの DCT-II に分解することで導かれる高速コサイン変換を 示します．この考えに基づく高速コサイン変換は B.G. Lee [参考文献] による 高速コサイン変換として知られているようです．
まず，DCT-II :
N-1
Ａ ＝ Σ ａ cos(π(j+1/2)k/N)
k j=0 j
N-1 (j+1/2)k
＝ Σ ａ Ｃ ， 0 <= k < N
j=0 j N
を次のように分解します．
N/2-1 (j+1/2)k
Ａ ＝ Σ （ａ ＋ ａ ） Ｃ
2k j=0 j N-1-j N/2
N/2-1 j+1/2 (j+1/2)k
Ａ ＋ Ａ ＝ Σ （ａ − ａ ）2 Ｃ Ｃ
2k-1 2k+1 j=0 j N-1-j N N/2
ここで，Ｃ は
j
Ｃ ＝ cos(πj/N)
N
です．この分解は，基数 2 の周波数間引き FFT の分解と似ていますが， Ａ の奇数の項の分解が異なり，簡単な漸化式になっています．したがって， 通常のバタフライ演算以外に後処理が必要になります．具体的には， まず分解の式に K=0 を代入し，1/2 倍して初項 Ａ_1 を求めます (Ａ_-1 = Ａ_1 という関係から)．次に，逐次的に Ａ の奇数項を求めて行きます． この後処理はスケーリングを除外すれば加減算だけで計算できます．しかし， これはバタフライの段数-1 回行わなければなりません．
次に，この方法による DCT-II のルーチンを示します．
| リスト2.3.3-3. B.G. Lee による 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[])
{
int m, mh, i, j, k, irev, j0, j1;
double theta, c2, xr;
/* ---- scrambler ---- */
i = 0;
for (j = 1; j < n - 1; j++) {
for (k = n >> 1; k > (i ^= k); k >>= 1);
if (j < i) {
xr = a[j];
a[j] = a[i];
a[i] = xr;
}
}
theta = 2 * atan(1.0) / n;
/* ---- butterflies ---- */
for (mh = 1; (m = mh << 1) <= n; mh = m) {
irev = 0;
for (i = 0; i < n; i += m) {
c2 = 2 * cos(theta * (irev + mh));
for (k = n >> 1; k > (irev ^= k); k >>= 1);
for (j = 0; j < mh; j++) {
j0 = i + j;
j1 = n - mh - i + j;
xr = a[j0] - a[j1];
a[j0] += a[j1];
a[j1] = c2 * xr;
}
}
}
/* ---- adds ---- */
for (m = n >> 1; (mh = m >> 1) >= 1; m = mh) {
for (j = 0; j < mh; j++) {
j0 = m + mh + j;
a[j0] = -a[j0] - a[j0 - m];
}
for (i = (m << 1) + mh; i < n; i += (m << 1)) {
for (j = 0; j < mh; j++) {
j0 = i + j;
j1 = j0 + m;
a[j0] -= a[j0 - m];
a[j1] = -a[j1] - a[j0];
}
}
}
/* ---- scaling ---- */
for (j = 1; j < n; j++) {
a[j] *= 0.5;
}
}
|
このルーチンの N=8 のときのデータフローは次のようになります．
 |
| 図2.3.3-2. B.G. Lee による FCT-II |
|---|
次に，演算量について少し考えます．1/2 によるスケーリングは， バタフライ演算の回転因子に含めることができます．したがって， 長さ N の DCT-II の計算に必要な乗算回数 μ と加算回数 α は次のように なります．
N
μ ＝ --- log N
2
3N
α ＝ ---- log N - N + 1
2
この演算量は非常に少なく，拡張 FFT の計算量を半分にする第一の方法に 対して，Split-Radix FFT を用い，さらに複素乗算を 3 回の実数乗算と 3 回の 実数加算で計算したときの演算量に匹敵します．また，この B.G. Lee の方法に 関しても大きい基数や Split-Radix にすれば演算量は減らせると思うかも しれませんが，実際はそう甘くはなく演算量はまったく同じになります． 基数をいじることは，ループ展開による高速化技法としてしか意味を持たなく なります．しかし，二次元 DCT を計算する場合は通常の FFT 同様に，後に示す Vector-Radix FFT を用いれば 乗算回数だけをさらに 3/4 に減らすことが可能になります．
DCT-III に関しても，この手続きを逆にすることで得られます． このとき，前処理は単なる和で直接計算できます．しかし，バタフライの 回転因子には 1/(2*cos(π*j/N)) と除算が入ってしまい，除算を 避けるようにあらかじめ計算しておかなければなりません．
最後に，この方法の欠点を挙げておきます．それは，バタフライ演算と 後処理の漸化式の計算で，長さ N に比例する丸め誤差が入るということです． したがって，あまり大きな計算には向きません．
ここでは，長さ N の DCT-II を一つの N/2 の長さの DCT-II と二つの N/4 の長さの DCT-II に分解することによる高速コサイン変換を示します．
まず，DCT-II :
N-1
Ａ ＝ Σ ａ cos(π(j+1/2)k/N) ， 0 <= k < N
k j=0 j
を次のように N/2 の長さの DCT-II と DCT-IV に分解します．
N/2-1
Ａ ＝ Σ （ａ ＋ ａ ） cos(π(j+1/2)k/(N/2))
2k j=0 j N-1-j
N/2-1
Ａ ＝ Σ （ａ − ａ ）cos(π(j+1/2)(k+1/2)/(N/2))
2k+1 j=0 j N-1-j
次に，N/2 の長さの DCT-IV を 2.3.2 で示した DCT-IV に対する分解を 用いて二つの N/4 の長さの DCT-II にします．この分解を繰り返すことで 高速算法が得られます．このとき同時に，DCT-IV に対する高速算法も 得られます．この高速算法は Split-Radix FFT とよく似た構造で L 字型の バタフライとなり，log(N) 段の演算は一段ずつ終了しなくなります． さらに，加減算からなる後処理が必要で多少複雑になります．また， この方法の演算量は DCT-IV の分解に現れる複素乗算を 3 回の加算と 3 回の乗算に置き換えたとき，B.G. Lee の方法とほぼ同じになります．
この方法の分解は，N 点の DCT を通常の DFT だと思って， Cooley-Tukey 型 FFT の分解を素直に用いたものです．この分解を表にすると 次のようになります．
| リスト2.3.3-4. DCT の分解 | |
|---|---|
| N+1 点 DCT-I | N+1 回の加算により N/2+1 点 DCT-I と N/2 点 DCT-III に分解 |
| N 点 DCT-II | N 回の加算により N/2 点 DCT-II と N/2 点 DCT-IV に分解 |
| N 点 DCT-IV | N-2 回の加算と N/2 回の複素数乗算により二組の N/2 点 DCT-II に分解 |
| N+1 点 DCT-I | N+1 回の加算により N/2+1 点 DCT-I と N/2 点 DCT-II に分解 |
| N 点 DCT-III | N 回の加算により N/2 点 DCT-III と N/2 点 DCT-IV に分解 |
| N 点 DCT-IV | N-2 回の加算と N/2 回の複素数乗算により二組の N/2 点 DCT-III に分解 |
この表から，長さ N の DCT を単純に半分の長さに分解したとき，別の型の 変換が混ざり，通常のDFT のように自分自身の変換には戻ってくれません． しかし，自分自身の変換に戻らなくても，この分解は繰り返すことができて， 四つのタイプのすべてが高速に計算できます．