はじめに
Fast cosine transformを説明します。B.G. Lee のPaperを参照します。
B. G. Lee A new algorithm to compute the discrete cosine transform
inverse DCT II
定義式は次の通りです。簡略化のため$C_{2N}^{(2n+1)k} = \cos \frac{\pi}{2N} (2n+1)k$と表記します。
\begin{eqnarray}
x_n &=& \sum_{k=0}^{N-1} X_k C_{2N}^{(2n+1)k} \\
\end{eqnarray}
kについて偶数と奇数に分けます。
\begin{eqnarray}
x_n &=&
\left\{
\begin{array}{ll}
g_n &=& \sum_{k=0}^{N/2-1} X_{2k} C_{2N}^{(2n+1)2k} \\
h'_n &=& \sum_{k=0}^{N/2-1} X_{2k+1} C_{2N}^{(2n+1)(2k+1)} \\
\end{array}
\right. \\
\end{eqnarray}
nについて前半 N/2 と 後半 N/2 に分けます。
次が成り立ちます。n=0 ... N/2-1 です。
\begin{eqnarray}
x_n &=& g_n + h'_n \\
x_{N-1-n} &=& g_n - h'_n \\
\end{eqnarray}
上記を示します。前半は自明です。後半 N/2 が成り立つことを説明します。
\begin{eqnarray}
x_{N-1-n} &=& \sum_{k=0}^{N/2-1} X_k C_{2N}^{(2(N-1-n)+1)k} \\
C_{2N}^{(2(N-1-n)+1)k} &=& C_{2N}^{2Nk-(2n+1)k} = \cos \left( k \pi - \frac{\pi (2n+1)k}{2N} \right) \\
\end{eqnarray}
次の関係式から後半が成り立つことが分かります。
\begin{eqnarray}
\cos \left(n \pi - \theta \right) =
\left\{
\begin{array}{ll}
+ \cos \theta & (n = even) \\
- \cos \theta & (n = odd) \\
\end{array}
\right. \\
\end{eqnarray}
次に再帰的に処理を行うために $g_n,h'_n$ を $x_n$と同じ式の形にします。
$x_n$は次の形でした。
\begin{eqnarray}
x_n &=& \sum_{k=0}^{N-1} X_{k} C_{2N}^{(2n+1)k} \\
\end{eqnarray}
$g_n$について考えます。
\begin{eqnarray}
g_n &=& \sum_{k=0}^{N/2-1} X_{2k} C_{2N}^{(2n+1)2k} \\
&=& \sum_{k=0}^{N/2-1} X_{2k} C_{2(N/2)}^{(2n+1)k} \\
G_k &=& X_{2k} \\
g_n &=& \sum_{k=0}^{N/2-1} G_k C_{2(N/2)}^{(2n+1)k} \\
\end{eqnarray}
$h'_n$について考えます。
\begin{eqnarray}
h'_n &=& \sum_{k=0}^{N/2-1} X_{2k+1} C_{2N}^{(2n+1)(2k+1)} \\
\end{eqnarray}
次の公式を使うことで、
\begin{eqnarray}
2\cos \alpha \cos \beta = \cos (\alpha + \beta) + \cos (\alpha - \beta)
\end{eqnarray}
次のように式変形できます。
\begin{eqnarray}
2 C_{2N}^{(2n+1)} C_{2N}^{(2n+1)(2k+1)} &=& C_{2N}^{(2n+1)2k} + C_{2N}^{(2n+1)2(k+1)} \\
\end{eqnarray}
上記を$h'_n$に適用します。
\begin{eqnarray}
2 C_{2N}^{(2n+1)} h'_n &=& \sum_{k=0}^{N/2-1} X_{2k+1} \left( C_{2N}^{(2n+1)2k} + C_{2N}^{(2n+1)2(k+1)} \right) \\
\end{eqnarray}
上記右辺2項を次のように変形します。
右肩の$2(k+1)$ を $2k$ にするために、2つ前の$X$を使います。
\begin{eqnarray}
\sum_{k=0}^{N/2-1} X_{2k+1} C_{2N}^{(2n+1)2(k+1)} &=& \sum_{k=0}^{N/2-1} X_{2k-1} C_{2N}^{(2n+1)2k} \\
\end{eqnarray}
\begin{eqnarray}
X_{-1} &=& 0 \\
C_{2N}^{(2n+1)2(N/2)} &=& C_{2}^{(2n+1)} = 0 \\
\end{eqnarray}
上記を$h'_n$に適用します。
\begin{eqnarray}
2 C_{2N}^{(2n+1)} h'_n &=& \sum_{k=0}^{N/2-1} \left( X_{2k+1} + X_{2k-1} \right) C_{2N}^{(2n+1)2k} \\
h'_n &=& \frac{1}{2 C_{2N}^{(2n+1)}} \sum_{k=0}^{N/2-1} H_k C_{2(N/2)}^{(2n+1)k} = \frac{1}{2 C_{2N}^{(2n+1)}} h_n \\
h_n &=& \sum_{k=0}^{N/2-1} H_k C_{2(N/2)}^{(2n+1)k} \\
H_n &=& X_{2k+1} + X_{2k-1} \\
\end{eqnarray}
式をまとめます。$x_n$は次の通りです。
\begin{eqnarray}
x_n &=& \sum_{k=0}^{N-1} X_k C_{2N}^{(2n+1)k} \\
n &=& 0,..., N-1 \\
\end{eqnarray}
$g_n,h_n$は次の通りです。
\begin{eqnarray}
g_n &=& \sum_{k=0}^{N/2-1} G_k C_{2(N/2)}^{(2n+1)k} \\
h_n &=& \sum_{k=0}^{N/2-1} H_k C_{2(N/2)}^{(2n+1)k} \\
n &=& 0,..., N/2-1 \\
\end{eqnarray}
$G_k,H_k$を$X_k$で表します。
\begin{eqnarray}
G_k &=& X_{2k} \\
H_k &=& X_{2k+1} + X_{2k-1} \\
\end{eqnarray}
$x_k$を$g_k,h_k$で表します。
\begin{eqnarray}
x_n &=& g_n + \frac{1}{2 C_{2N}^{(2n+1)}} h_n \\
x_{N-1-n} &=& g_n - \frac{1}{2 C_{2N}^{(2n+1)}} h_n \\
\end{eqnarray}
次にフロー図を示します。N=8の場合です。
forward DCT II
定義式は次の通りです。
\begin{eqnarray}
X_k &=& \sum_{n=0}^{N-1} x_n C_{2N}^{(2n+1)k} \\
\end{eqnarray}
nについて前半 N/2 と 後半 N/2 に分けます。
次が成り立ちます。n=0 ... N/2-1 です。
\begin{eqnarray}
X_k &=& \sum_{n=0}^{N/2-1} x_n C_{2N}^{(2n+1)k} + \sum_{n=0}^{N/2-1} x_{N-1-n} C_{2N}^{(2(N-1-n)+1)k} \\
\end{eqnarray}
kについて偶数と奇数に分けます。
\begin{eqnarray}
X_{2k}
&=& \sum_{n=0}^{N/2-1} x_n C_{2N}^{(2n+1)2k} + \sum_{n=0}^{N/2-1} x_{N-1-n} C_{2N}^{(2(N-1-n)+1)2k} \\
&=& \sum_{n=0}^{N/2-1} x_n C_{2N/2}^{(2n+1)k} + \sum_{n=0}^{N/2-1} x_{N-1-n} C_{2N/2}^{(2n+1)k} \\
&=& \sum_{n=0}^{N/2-1} (x_n - x_{N-1-n}) C_{2N/2}^{(2n+1)k} = \sum_{n=0}^{N/2-1} g_n C_{2N/2}^{(2n+1)k} = G_k \\
X_{2k+1}
&=& \sum_{n=0}^{N/2-1} x_n C_{2N}^{(2n+1)(2k+1)} + \sum_{n=0}^{N/2-1} x_{N-1-n} C_{2N}^{(2(N-1-n)+1)(2k+1)} \\
&=& \sum_{n=0}^{N/2-1} x_n C_{2N}^{(2n+1)(2k+1)} - \sum_{n=0}^{N/2-1} x_{N-1-n} C_{2N}^{(2n+1)(2k+1)} \\
&=& \sum_{n=0}^{N/2-1} (x_n - x_{N-1-n}) C_{2N}^{(2n+1)(2k+1)} \\
&=& \sum_{n=0}^{N/2-1} \frac{(x_n - x_{N-1-n})}{2 C_{2N}^{(2n+1)}} \left (C_{2N/2}^{(2n+1)(k+1)} + C_{2N/2}^{(2n+1)k} \right) \\
&=& \sum_{n=0}^{N/2-1} h_n \left (C_{2N/2}^{(2n+1)(k+1)} + C_{2N/2}^{(2n+1)k} \right) = H_{k+1} + H_k \\
\end{eqnarray}
式をまとめます。$X_n$は次の通りです。
\begin{eqnarray}
X_k &=& \sum_{n=0}^{N-1} x_n C_{2N}^{(2n+1)k} \\
k &=& 0, ..., N-1 \\
\end{eqnarray}
$G_n, H_n$は次の通りです。
\begin{eqnarray}
G_k &=& \sum_{n=0}^{N/2-1} g_n C_{2N/2}^{(2n+1)k} \\
H_k &=& \sum_{n=0}^{N/2-1} h_n C_{2N/2}^{(2n+1)k} \\
k &=& 0, ..., N/2-1 \\
\end{eqnarray}
$g_n, h_n$を$x_n$で表します。
\begin{eqnarray}
g_n &=& x_n + x_{N-1-n} \\
h_n &=& \frac{x_n - x_{N-1-n}}{2 C^{2n+1}_{2N}} \\
\end{eqnarray}
$X_k$を$G_k,H_k$で表します。
\begin{eqnarray}
X_{2k} &=& G_k \\
X_{2k+1} &=& H_k + H_{k+1} \\
\end{eqnarray}
次にフロー図を示します。N=8の場合です。
DCT IV
DCT IVはDCT IIに変換できます。次のようにして変換します。
\begin{eqnarray}
x_n
&=& \sum_{k=0}^{N-1}{X_k C_{4N}^{(2n+1)(2k+1)}} \\
&=& \frac{1}{2 C_{4N}^{2n+1}} \sum_{k=0}^{N-1}{X_k \left( C_{2N}^{(2n+1)k} + C_{2N}^{(2n+1)(k+1)} \right) } \\
&=& \frac{1}{2 C_{4N}^{2n+1}} \sum_{k=0}^{N-1}{(X_{k-1} + X_k) C_{2N}^{(2n+1)k}} \\
&=& \frac{1}{2 C_{4N}^{2n+1}} y_n \\
y_n &=& \sum_{k=0}^{N-1}{Y_k C_{2N}^{(2n+1)k}} \\
Y_k &=& X_{k-1} + X_k \\
\end{eqnarray}
code example
makefile
CFLAGS=-I. -Wall -Werror -O2 -march=native
INCS=
OBJS=test.o
LIBS=
TARGET=test
%.o: %.c $(INCS)
$(CC) $(CFLAGS) -c -o $@ $<
$(TARGET): $(OBJS)
$(CC) $(CFLAGS) -o $@ $^ $(LIBS)
clean:
rm -rf $(TARGET) *.o
test.c
# include <math.h>
# include <stdio.h>
# include <string.h>
const double C4[] = {
7.07106781186548E-01,
};
const double C8[] = {
5.41196100146197E-01, // 1/(2*cos(pi/8*1))
1.30656296487638E+00, // 1/(2*cos(pi/8*3))
};
const double C16[] = {
5.09795579104159E-01, // 1/(2*cos(pi/16*1))
6.01344886935045E-01, // 1/(2*cos(pi/16*3))
8.99976223136416E-01, // 1/(2*cos(pi/16*5))
2.56291544774151E+00, // 1/(2*cos(pi/16*7))
};
const double C32[] = {
5.02419286188156E-01, // 1/(2*cos(pi/32*1))
5.22498614939689E-01, // 1/(2*cos(pi/32*3))
5.66944034816358E-01, // 1/(2*cos(pi/32*5))
6.46821783359990E-01, // 1/(2*cos(pi/32*7))
7.88154623451250E-01, // 1/(2*cos(pi/32*9))
1.06067768599035E+00, // 1/(2*cos(pi/32*11))
1.72244709823833E+00, // 1/(2*cos(pi/32*13))
5.10114861868916E+00, // 1/(2*cos(pi/32*15))
};
static void idctII2(double* input, double* output)
{
double* x = output;
double* X = input;
double tmp = C4[0] * X[1];
x[0] = X[0] + tmp;
x[1] = X[0] - tmp;
}
static void idctII4(double* input, double* output)
{
double* x = output;
static double in[2], g[2], h[2], tmp[2];
in[0] = input[0];
in[1] = input[2];
idctII2(in, g);
in[0] = input[1];
in[1] = input[3] + input[1];
idctII2(in, h);
tmp[0] = C8[0] * h[0];
tmp[1] = C8[1] * h[1];
x[0] = g[0] + tmp[0];
x[1] = g[1] + tmp[1];
x[2] = g[1] - tmp[1];
x[3] = g[0] - tmp[0];
}
static void idctII8(double* input, double* output)
{
double* x = output;
static double in[4], g[4], h[4], tmp[4];
in[0] = input[0];
in[1] = input[2];
in[2] = input[4];
in[3] = input[6];
idctII4(in, g);
in[0] = input[1];
in[1] = input[3] + input[1];
in[2] = input[5] + input[3];
in[3] = input[7] + input[5];
idctII4(in, h);
tmp[0] = C16[0] * h[0];
tmp[1] = C16[1] * h[1];
tmp[2] = C16[2] * h[2];
tmp[3] = C16[3] * h[3];
x[0] = g[0] + tmp[0];
x[1] = g[1] + tmp[1];
x[2] = g[2] + tmp[2];
x[3] = g[3] + tmp[3];
x[4] = g[3] - tmp[3];
x[5] = g[2] - tmp[2];
x[6] = g[1] - tmp[1];
x[7] = g[0] - tmp[0];
}
static void fdctII2(double* input, double* output)
{
double* X = output;
double* x = input;
X[0] = x[0] + x[1];
X[1] = (x[0] - x[1]) * C4[0];
}
static void fdctII4(double* input, double* output)
{
double* X = output;
static double in[2], G[2], H[2];
in[0] = input[0] + input[3];
in[1] = input[1] + input[2];
fdctII2(in, G);
in[0] = C8[0] * (input[0] - input[3]);
in[1] = C8[1] * (input[1] - input[2]);
fdctII2(in, H);
X[0] = G[0];
X[1] = H[0] + H[1];
X[2] = G[1];
X[3] = H[1];
}
static void fdctII8(double* input, double* output)
{
double* X = output;
static double in[4], G[4], H[4];
in[0] = input[0] + input[7];
in[1] = input[1] + input[6];
in[2] = input[2] + input[5];
in[3] = input[3] + input[4];
fdctII4(in, G);
in[0] = C16[0] * (input[0] - input[7]);
in[1] = C16[1] * (input[1] - input[6]);
in[2] = C16[2] * (input[2] - input[5]);
in[3] = C16[3] * (input[3] - input[4]);
fdctII4(in, H);
X[0] = G[0];
X[1] = H[0] + H[1];
X[2] = G[1];
X[3] = H[1] + H[2];
X[4] = G[2];
X[5] = H[2] + H[3];
X[6] = G[3];
X[7] = H[3];
}
static void dctIV8(double* input, double* output)
{
int i;
static double in[8];
in[0] = input[0];
for (i=1; i<8; i++) {
in[i] = input[i-1] + input[i];
}
idctII8(in, output);
for (i=0; i<8; i++) {
output[i] *= C32[i];
}
}
static void print_data(const char* str, const double *data, int N)
{
int i;
printf("%s", str);
for (i=0; i<N; i++) {
printf("%+1.3e ", data[i]);
}
printf("\n");
}
int main(int argc, char* argv[])
{
double input[] = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0};
double output[8];
print_data("input\n", input, 8);
memset(output, 0, sizeof(output));
fdctII8(input, output);
print_data("fdctII8\n", output, 8);
memset(output, 0, sizeof(output));
idctII8(input, output);
print_data("idctII8\n", output, 8);
memset(output, 0, sizeof(output));
dctIV8(input, output);
print_data("dctIV8\n", output, 8);
return 0;
}
console
$ make;./test.exe
cc -I. -Wall -Werror -O2 -march=native -c -o test.o test.c
cc -I. -Wall -Werror -O2 -march=native -o test test.o
input
+0.000e+000 +1.000e+000 +2.000e+000 +3.000e+000 +4.000e+000 +5.000e+000 +6.000e+000 +7.000e+000
fdctII8
+2.800e+001 -1.288e+001 +0.000e+000 -1.347e+000 +0.000e+000 -4.018e-001 +0.000e+000 -1.014e-001
idctII8
+1.459e+001 -1.615e+001 +6.358e+000 -5.495e+000 +2.864e+000 -2.459e+000 +9.404e-001 -6.464e-001
dctIV8
+1.236e+001 -1.576e+001 +6.963e+000 -6.391e+000 +4.586e+000 -4.404e+000 +3.839e+000 -3.793e+000