Discrete cosine transformを理解する

はじめに

Discrete cosine transform(DCT)について説明します。
DCTは映像、音声圧縮に利用されます。
いくつかのバリエーションがありますが、ここではDCT II /DCT IVについて説明します。

DCT II

forward

X_k = \sqrt{\frac{2}{N}} C_k \sum_{n=0}^{N-1}{x_n \cos \left[ \frac{\pi}{2N} (2n+1)k \right] }

inverse

x_n = \sqrt{\frac{2}{N}} \sum_{k=0}^{N-1}{C_k X_k \cos \left[ \frac{\pi}{2N} (2n+1)k \right] }

$C_k$はk=0 で $\frac{1}{\sqrt{2}}$、k!=0 で 1となります。

DCT IV

forward

X_k = \sqrt{\frac{2}{N}} \sum_{n=0}^{N-1}{x_n \cos \left[ \frac{\pi}{4N} (2n+1)(2k+1) \right] }

inverse

x_n = \sqrt{\frac{2}{N}} \sum_{k=0}^{N-1}{X_k \cos \left[ \frac{\pi}{4N} (2n+1)(2k+1) \right] }

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>

static void fdctII(const double *input, double *output, int N);
static void idctII(const double *input, double *output, int N);
static void dctIV(const double *input, double *output, int N);
static void print_data(const char* str, const double *data, int N);

int main(int argc, char* argv[])
{
  double input[] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
  double output1[8], output2[8];

  print_data("input\n", input, 8);

  memset(output1, 0, sizeof(output1));
  fdctII(input,   output1, 8);
  print_data("fdctII\n", output1, 8);
  memset(output2, 0, sizeof(output2));
  idctII(output1, output2, 8);
  print_data("idctII\n", output2, 8);

  memset(output1, 0, sizeof(output1));
  dctIV(input,   output1, 8);
  print_data("fdctIV\n", output1, 8);
  memset(output2, 0, sizeof(output2));
  dctIV(output1, output2, 8);
  print_data("idctIV\n", output2, 8);

  return 0;
}

static void fdctII(const double *input, double *output, int N)
{
  int n, k;
  double c1 = sqrt(2.0/(double)N);
  double c2 = 1.0 / sqrt(2.0);

  for (k=0; k<N; k++) {
    output[k] = 0.0;
    for (n=0; n<N; n++) {
      output[k] += input[n] * cos(M_PI/2.0/(double)N*(2.0*n+1.0)*k);
    }
    output[k] *= c1;
    if (k == 0) {
      output[k] *= c2;
    }
  }
}

static void idctII(const double *input, double *output, int N)
{
  int n, k;
  double c1 = sqrt(2.0/(double)N);
  double c2 = 1.0 / sqrt(2.0);

  for (n=0; n<N; n++) {
    output[n] = 0.0;
    for (k=0; k<N; k++) {
      if (k == 0) {
        output[n] += c2 * input[k] * cos(M_PI/2.0/(double)N*(2.0*(double)n+1.0)*(double)k);
      } else {
        output[n] +=      input[k] * cos(M_PI/2.0/(double)N*(2.0*(double)n+1.0)*(double)k);
      }
    }
    output[n] *= c1;
  }
}

static void dctIV(const double *input, double *output, int N)
{
  int n, k;
  double c1 = sqrt(2.0/(double)N);

  for (n=0; n<N; n++) {
    output[n] = 0.0;
    for (k=0; k<N; k++) {
      output[n] += input[k] * cos(M_PI/4.0/(double)N*(2.0*(double)n+1.0)*(2.0*(double)k+1.0));
    }
    output[n] *= c1;
  }
}

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");
}

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
+1.000e+000 +2.000e+000 +3.000e+000 +4.000e+000 +5.000e+000 +6.000e+000 +7.000e+000 +8.000e+000
fdctII
+1.273e+001 -6.442e+000 -1.332e-015 -6.735e-001 +0.000e+000 -2.009e-001 -1.577e-014 -5.070e-002
idctII
+1.000e+000 +2.000e+000 +3.000e+000 +4.000e+000 +5.000e+000 +6.000e+000 +7.000e+000 +8.000e+000
fdctIV
+8.732e+000 -8.740e+000 +4.012e+000 -3.590e+000 +2.616e+000 -2.485e+000 +2.181e+000 -2.148e+000
idctIV
+1.000e+000 +2.000e+000 +3.000e+000 +4.000e+000 +5.000e+000 +6.000e+000 +7.000e+000 +8.000e+000

Appendix

inverse DCT II

inverse DCT II により$x_n$が復元できることを確認します。
inverse DCT IIの式にforward DCT IIを代入します。

\begin{eqnarray}
&& \sqrt{\frac{2}{N}} \sum_{k=0}^{N-1}{C_k X_k \cos \left[ \frac{\pi}{2N} (2n+1)k \right] } \\
&=& \sqrt{\frac{2}{N}} \sum_{k=0}^{N-1}{C_k \left [ \sqrt{\frac{2}{N}} C_k \sum_{m=0}^{N-1}{x_m \cos \left[ \frac{\pi}{2N} (2m+1)k \right] } \right] \cos \left[ \frac{\pi}{2N} (2n+1)k \right] } \\
&=& \frac{2}{N} \sum_{k=0}^{N-1}{C_k^2 \sum_{m=0}^{N-1}{x_m \cos \left[ \frac{\pi}{2N} (2m+1)k \right] } \cos \left[ \frac{\pi}{2N} (2n+1)k \right] } \\
\end{eqnarray} \\

$\sum$の順番を入れ替えてまとめます。

\frac{2}{N} \sum_{m=0}^{N-1}{x_m \sum_{k=0}^{N-1}{C_k^2 \cos \left[ \frac{\pi}{2N} (2m+1)k \right] } \cos \left[ \frac{\pi}{2N} (2n+1)k \right] }

次の公式を使って

\cos \alpha \cos \beta = \frac{1}{2} \left( \cos(\alpha + \beta) + \cos(\alpha - \beta) \right)

次のように式変形します。

\begin{eqnarray} \\
\cos \left[ \frac{\pi}{2N} (2m+1)k \right] \cos \left[ \frac{\pi}{2N} (2n+1)k \right] &=& \frac{1}{2} \left[ \cos \left[ \frac{\pi}{N} (m+n+1)k \right] + \cos \left[ \frac{\pi}{N} (m-n)k \right] \right]
\end{eqnarray} \\

上記を使って次のように式変形します。式1とします。

\frac{1}{N} \sum_{m=0}^{N-1}{x_m \sum_{k=0}^{N-1}{C_k^2 } \left[ \cos \left[ \frac{\pi}{N} (m+n+1)k \right] + \cos \left[ \frac{\pi}{N} (m-n)k \right] \right] }

次の公式があります。公式1とします。

\sum_{k=0}^{N-1}{\cos (k \theta + \phi )} = \frac{\sin (\frac{N \theta}{2}) \cos(\frac{(N-1)\theta}{2} + \phi )}{\sin \frac{\theta}{2}}

上記公式を使い$\theta = \frac{\pi}{N} l, \phi = 0$ として次の式を考えます。式2とします。

\sum_{k=0}^{N-1}{\cos (\frac{\pi}{N} l k)} = \frac{\sin (\frac{\pi}{2} l) \cos(\frac{\pi}{2} l - \frac{\pi}{2N} l)}{\sin \frac{\pi}{2N} l}

上記式は$l$が偶数であれば、$\sin \frac{\pi}{2} l = 0$のため0になります。
奇数であれば1になります。
$l=1$では次のようになります。

\frac{\sin (\frac{\pi}{2}) \cos(\frac{\pi}{2} - \frac{\pi}{2N})}{\sin \frac{\pi}{2N}} = \frac{1 \sin \frac{\pi}{2N}}{\sin \frac{\pi}{2N}} = 1

$l=3$では次のようになります。

\frac{\sin (\frac{3 \pi}{2}) \cos(\frac{3 \pi}{2} - \frac{3 \pi}{2N})}{\sin \frac{3 \pi}{2N}} = \frac{(-1)(- \sin \frac{3 \pi}{2N})}{\sin \frac{3 \pi}{2N}} = 1

式1の次の部分について考えます。式3とします。

\sum_{k=0}^{N-1} \left[ \cos \left[ \frac{\pi}{N} (m+n+1)k \right] + \cos \left[ \frac{\pi}{N} (m-n)k \right] \right]

式2が2つ並んだ形をしています。$l=m+n+1 $ or $m-n$ となります。
m+n+1が偶数であればm-nは奇数になります。
逆にm+n+1が奇数であればm-nは偶数になります。
また、式3の2項目はm=nの時にNになります。m!=nでは0になります。
式3の1項はm=nの時、m+n+1は奇数となるため、1となります。

したがって式3は次のようになります。$\delta_{mn}$はクロネッカーのデルタです。m=nで1,m!=nで0となります。

\sum_{k=0}^{N-1} \left[ \cos \left[ \frac{\pi}{N} (m+n+1)k \right] + \cos \left[ \frac{\pi}{N} (m-n)k \right] \right] = 1 + N \delta_{mn}

$C_k^2$について考えます。k=0の時に1/2になります。
次のように$C_k^2$がある場合はk=0 で 1となります。
ない場合は2のため、これは式3を-1していることになります。

\frac{1}{2} \left[ \cos \left[ \frac{\pi}{N}(m+n+1) 0) \right] + \cos \left[ \frac{\pi}{N}(m-n) 0) \right] \right] = \frac{1}{2} \left[ \cos 0 + \cos 0  \right] = 1

したがって次のようになります。

\sum_{k=0}^{N-1} C_{k}^{2} \left[ \cos \left[ \frac{\pi}{N} (m+n+1)k \right] + \cos \left[ \frac{\pi}{N} (m-n)k \right] \right] = N \delta_{mn}

上記を式1に適用します。

\begin{eqnarray}
\frac{1}{N} \sum_{m=0}^{N-1}{x_m \sum_{k=0}^{N-1}{C_k^2 } \left[ \cos \left[ \frac{\pi}{N} (m+n+1)k \right] + \cos \left[ \frac{\pi}{N} (m-n)k \right] \right] } = \frac{1}{N} \sum_{m=0}^{N-1}{x_m} N \delta_{mn} = x_n
\end{eqnarray}

よって、inverse DCT II により$x_n$が復元できることがわかりました。

inverse DCT IV

inverse DCT IV により$x_n$が復元できることを確認します。
考え方はinverse DCT IIと同じです。

inverse DCT IVの式にforward DCT IVを代入します。

\begin{eqnarray}
&& \sqrt{\frac{2}{N}} \sum_{k=0}^{N-1}{X_k \cos \left[ \frac{\pi}{4N} (2n+1)(2k+1) \right] } \\
&=& \sqrt{\frac{2}{N}} \sum_{k=0}^{N-1}{\left [ \sqrt{\frac{2}{N}} \sum_{m=0}^{N-1}{x_m \cos \left[ \frac{\pi}{4N} (2m+1)(2k+1) \right] } \right] \cos \left[ \frac{\pi}{4N} (2n+1)(2k+1) \right] } \\
&=& \frac{2}{N} \sum_{k=0}^{N-1}{\sum_{m=0}^{N-1}{x_m \cos \left[ \frac{\pi}{4N} (2m+1)(2k+1) \right] } \cos \left[ \frac{\pi}{4N} (2n+1)(2k+1) \right] } \\
&=& \frac{2}{N} \sum_{m=0}^{N-1}x_m{\sum_{k=0}^{N-1}{\cos \left[ \frac{\pi}{4N} (2m+1)(2k+1) \right] } \cos \left[ \frac{\pi}{4N} (2n+1)(2k+1) \right] } \\
&=& \frac{1}{N} \sum_{m=0}^{N-1}x_m{\sum_{k=0}^{N-1}\left[{ \cos \left[ \frac{\pi}{2N} (m+n+1)(2k+1) \right] } + \cos \left[ \frac{\pi}{2N} (m-n)(2k+1) \right] \right] } \\
\end{eqnarray}

上記最後の式を式4とします。
公式1を使い、$\theta = \frac{\pi}{N} l, \phi = \frac{\pi}{2N} l$ として次の式を考えます。式5とします。

\begin{eqnarray}
\sum_{k=0}^{N-1}{\cos (\frac{\pi}{N} l k + \frac{\pi}{2N}l)} &=& \frac{\sin (\frac{\pi}{2} l) \cos(\frac{\pi}{2N}l + \frac{\pi}{2} l - \frac{\pi}{2N} l)}{\sin \frac{\pi}{2N} l} = \frac{\sin (\frac{\pi}{2} l) \cos(\frac{\pi}{2} l)}{\sin \frac{\pi}{2N} l} \\
&=& \frac{1}{2} \frac{\sin (\pi l)}{\sin \frac{\pi}{2N}l}
\end{eqnarray}

$\sin (\pi l) = 0$のため、上記は0になります。

式5を式4に適用します。$l = $ m+n+1 or m-nとします。

\begin{eqnarray}
&& \frac{1}{N} \sum_{m=0}^{N-1}x_m{\sum_{k=0}^{N-1}\left[{ \cos \left[ \frac{\pi}{2N} (m+n+1)(2k+1) \right] } + \cos \left[ \frac{\pi}{2N} (m-n)(2k+1) \right] \right] } \\
&=& \frac{1}{N} \sum_{m=0}^{N-1} x_m N \delta_{mn} = x_n \\
\end{eqnarray}

よって、inverse DCT IV により$x_n$が復元できることがわかりました。