ここでの拡張 DFT は，時間方向または周波数方向の点をずらした 変換のことで，
N-1 (j+δ1)(k+δ2)
Ａ ＝ Σ ａ Ｗ
k j=0 j N
で定義します．この逆変換も拡張 DFT で
1 N-1 -(j+δ2)(k+δ1)
ａ ＝ --- Σ Ａ Ｗ
k N j=0 j N
となります．もし，δ1, δ2 が有理数ならば，拡張 DFT は通常の DFT を 適当な間隔で点をとばして計算したものに相当します．たとえば，δ1=1/2, δ2=0 の場合は長さ 2N の奇数点上での DFT になります．
この拡張 DFT は，通常の FFT の変換前のデータに Ｗ^δ2 j を掛けて， 変換後のデータに Ｗ^δ1 k を掛ければ計算できることになります． しかし，δ1 または δ2 のどちらかがゼロのこの変換は，通常の FFT の 三角関数データをずらすだけで容易に実行できます [参考文献]．したがって，これらの乗算の 片方は簡単に省略することができます．
例えば，時間方向に δ1 シフトした拡張 DFT は周波数間引き FFT の 三角関数を δ1 だけずらすことで計算できます．具体的なルーチンは 次のようになります．
| リスト2.2-1. 周波数間引きによる拡張 FFT |
|---|
/*
F[k]=Σ_j=0^n-1 a[j]*exp(±2*pi*i*(j+delta1)*k/n),0<=k<n を計算する.
n はデータ数で 2 の整数乗, theta は ±2*PI/n,
ar[0...n-1], ai[0...n-1], はデータの実部, 虚部.
*/
#include <math.h>
void gfft1(int n, double theta, double delta1, double ar[], double ai[])
{
int m, mh, i, j, k;
double wr, wi, xr, xi;
for (m = n; (mh = m >> 1) >= 1; m = mh) {
for (i = 0; i < mh; i++) {
wr = cos(theta * (i + delta1));
wi = sin(theta * (i + delta1));
for (j = i; j < n; j += m) {
k = j + mh;
xr = ar[j] - ar[k];
xi = ai[j] - ai[k];
ar[j] += ar[k];
ai[j] += ai[k];
ar[k] = wr * xr - wi * xi;
ai[k] = wr * xi + wi * xr;
}
}
theta *= 2;
}
/* ---- unscrambler ---- */
i = 0;
for (j = 1; j < n - 1; j++) {
for (k = n >> 1; k > (i ^= k); k >>= 1);
if (j < i) {
xr = ar[j];
xi = ai[j];
ar[j] = ar[i];
ai[j] = ai[i];
ar[i] = xr;
ai[i] = xi;
}
}
}
|
同様に，周波数方向に δ2 だけシフトした拡張 DFT は時間間引き FFT の 三角関数を δ2 だけずらすことで計算できます．これらの拡張 DFT に対する FFT を，拡張 FFT と呼ぶことにします．
次に実数データに対する拡張 FFT の計算について 考えます．これは， 実数専用の FFT ルーチンを用いる方法の第二の方法 を用いて計算できます．例えば，実数データの時間 1/2 シフト拡張 FFT は 三角関数データを 1/2 だけずらすことで計算できます．このルーチンは 次のようになります．
| リスト2.2-2. 実数データの時間 1/2 シフト拡張 FFT |
|---|
/*
実数データの時間 1/2 シフト拡張 DFT :
F[k]=Σ_j=0^n-1 a[j]*exp(-2*pi*i*(j+1/2)*k/n),0<=k<n を計算する.
出力 a[0...n/2-1], a[n/2...n-1] はそれぞれ F[0...n/2-1] の実部,
F[n/2...n-1] の符号を付けかえた虚部に相当する．
*/
#include <math.h>
void grfft1h(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 = -8 * atan(1.0) / n;
for (mh = 1; (m = mh << 1) <= n; mh = m) {
/* ---- real to real butterflies ---- */
for (jr = 0; jr < n; jr += m) {
kr = jr + mh;
xr = a[jr] - a[kr];
a[jr] += a[kr];
a[kr] = xr;
}
if (mh == 1) continue;
/* ---- real to complex butterflies ---- */
mq = mh >> 1;
irev = mq;
for (ji = mq; ji < n; ji += m) {
wr = cos(0.5 * theta * irev);
wi = sin(0.5 * theta * irev);
for (k = n >> 2; k > (irev ^= k); k >>= 1);
ki = ji + mh;
xr = a[ji];
xi = a[ki];
a[ji] = wr * xr + wi * xi;
a[ki] = wr * xi - wi * xr;
}
if (mh == 2) continue;
/* ---- complex to complex butterflies ---- */
irev = mq;
for (i = 0; i < n; i += m) {
wr = cos(theta * irev);
wi = sin(theta * irev);
for (k = n >> 2; k > (irev ^= k); k >>= 1);
for (j = 1; j < mq; j++) {
jr = i + j;
kr = jr + mh;
ki = i + mh - j;
ji = ki + mh;
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;
}
}
}
}
|
この実数データの時間 1/2 シフト拡張 FFT は，離散コサイン変換 DCT-II と 離散サイン変換 DST-II を同時に計算していることになります．したがって， このルーチンを二つに分離すれば，DCT-II と DST-II のルーチンが 作成できることになります．これは 次の 2.3.3で示します．