LoginSignup
2
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 5 years have passed since last update.

@tobira-code

Discrete cosine transformを理解する

Last updated at Posted at 2019-05-27

#はじめに
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$が復元できることがわかりました。

2
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?