研究室内向けに公開していたドキュメントです。

離散コサイン変換 DCTについてまとめてみた

本などで調べようにも、DCTはjpegの変換で用いられることから2次元の文献くらいしか残っていない。
DCT の起源である参考文献[1]をベースに、
改めて１次元のDCT変換方法について調べてみました。

離散コサイン変換 (DCT: Discrete Cosine Transform)

離散コサイン変換は[1]より画像処理に関連する問題である、パターン認識とWenerフィルタリングで活用するために開発されたらしいです。
今日では、jpegフォーマットへ圧縮する際に空間領域から周波数領域へ変換する時に用いられています。

そんなDCT でやっていることは、信号をcos波形の周期別に分解する作業をやっています(自己解釈)。
その結果どの周期のcos波形をどれだけ使っているのかを表しているのがDCT係数になります。
元の信号に戻す際には、DCT 係数のレベルに応じてcos信号を足し合わせれば元の信号に戻すことができます→IDCT。

DCT（時間領域→周波数領域）

$$X(0) =\frac{\sqrt{2}}{N} \sum\limits_{n=1}^{N-1} x(n) \tag{1}$$$$X(k) =\frac{2}{N} \sum\limits_{n=1}^{N-1} x(n)\cos{\frac{(2n+1)k\pi}{2N}}, k=1,2,...,(N-1) \tag{2}$$

DCTは、上記の式で周波数領域に変換することができます。
$x(n)$は、信号サンプル、$X(k)$は、DCT係数になります。
$N$は、サンプル数を表しており、DCT係数も同じの数になります。

$X(0)$は、直流成分になるので平均値を求めるような式になっています。
$X(0)$を除く$X(k)$は、$k$で示すcos周期信号が信号$x(n)$で使われているのかを算出する式になっています。

DCTでは上記2つの指揮を用いて時間領域から周波数領域へ変換を行っています。

式(1)の計算

式(1)直流成分に関しては、基本的に平均を求めれば良いのでcosはつけないで計算します。
プログラムを書く場合は以下のように書きます。

c dct1.c

#include <stdio.h>
#include <math.h>
int main(void) {
    int i, j; // loop
    double x[10] = {10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 15.0, 14.0, 13.0}; // 入力信号
    double coefficients[10]; // DCT係数
    const int size = 10; // 配列サイズ

    // 直流成分の計算
    coefficients[0] = 0;
    for (i = 0; i < size; i++) {
        coefficients[0] += x[i];
    }
    coefficients[0] *= sqrt(2) / size;
    
    return 0;
}

平均のプログラムにsqrt(2)追加するだけなので、比較的簡単に書けると思います。

式(2)の計算

次に式(2)の計算をしていきましょう。
プログラムを書く場合は以下のように書きます。

c dct2.c

#include <stdio.h>
#include <math.h>
int main(void) {
    int i, j; // loop
    double x[10] = {10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 15.0, 14.0, 13.0}; // 入力信号
    double coefficients[10]; // DCT係数
    const int size = 10; // 配列サイズ

    // 直流成分の計算
    { /* 省略 */ }

    // 周期成分の計算
    for (i = 1; i < size; i++) {
        coefficients[i] = 0;
        for (j = 0; j < size; j++) {
            double vector = (2 * j + 1) * i * M_PI / (2 * size); // cos内をあらかじめ計算
            coefficients[i] += x[j] * cos(vector);
        }
        coefficients[i] *= 2.0 / size;
    }
    
    return 0;
}

cos内をあらかじめ計算しておくのがおすすめです。
一気に計算する方が良さそうな気もしますが、分けて計算する方がプログラムの読みやすさは格段に上がります。

計算結果のアウトプット

計算結果を確認してみましょう。
以下は、結果を表示するためのプログラム例です。

c dct2.c

#include <stdio.h>
#include <math.h>
int main(void) {
    int i, j; // loop
    double x[10] = {10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 15.0, 14.0, 13.0}; // 入力信号
    double coefficients[10]; // DCT係数
    const int size = 10; // 配列サイズ

    // 直流成分の計算
    { /* 省略 */ }

    // 周期成分の計算
    { /* 省略 */ }

    // 計算結果の表示
    printf("signal -> coefficients\n");
    for (i = 0; i < size; i++) {
        printf("%6f,\t%6f\n", x[i], coefficients[i]);
    }
    
    return 0;
}

結果はこんな感じで表示されると思います。

signal -> coefficients
10.000000,	18.809040
11.000000,	-1.855162
12.000000,	-1.611496
13.000000,	0.479211
14.000000,	-0.323607
15.000000,	-0.141421
16.000000,	0.055503
15.000000,	-0.019705
14.000000,	-0.123607
13.000000,	0.091336
Program ended with exit code: 0

IDCT（周波数領域→時間領域）

次は、IDCTをやってみましょう。
$$x(n) =\frac{1}{\sqrt{2}}X(0)+ \sum\limits_{k=1}^{N-1} X(k)\cos{\frac{(2n+1)k\pi}{2N}}, n=0,1,...,(N-1) \tag{3}$$
上記の式はIDCTの式になります。
DCTとは違い、式は一つになっています。
式(1)がなくなったように見えるのですが、DCTの式(1)で計算した分はどこに行ったのかというと、第１項の部分
$\frac{1}{\sqrt{2}}X(0)$
がその式に当たります。
式(1)は直流成分を計算していましたので、式(1)で$\sqrt{2}$かけていた分、割り戻すことで元の値に戻るようにしています。
疑問が解けたところで、実際にプログラムを作っていきましょう。

式(3)の計算

式(3)を見ると基本的には足し合わせるだけで計算できそうです。
それぞれの項ごとに足し合わせるような形でプログラムを組んでいきましょう。
式(3)のプログラムをかくと以下のようになります。

c idct1.c

#include <stdio.h>
#include <math.h>
int main(void) {
    int i, j; // loop
    double coefficients[10] = {7.778175, -4.036036, -0.000000, -0.432302, -0.000000, -0.141421, -0.000000, -0.057185, 0.000000, -0.016036}; // DCT係数
    double resample[10]; // IDCTから生成される信号
    const int size = 10; // 配列サイズ
    
    // idct
    for (i = 0; i < size; i++) {
        resample[i] = coefficients[0] / sqrt(2); // 直流成分の計算
        for (j = 1; j < size; j++) {
            double vector = (2 * i + 1) * j * M_PI / (2 * size);
            resample[i] += coefficients[j] * cos(vector);
        }
    }
    
    return 0;
}

こんな感じに書いてみました。
上のように書けましたか？

DCTの時と同様、cosの中身を分けて書くと可読性の高いプログラムを書くことができます。
また、cosの中身であるvectorは、for文の変数に気をつければ、DCTと全く同じなので間違いにくくなります。

計算結果のアウトプット

計算結果を確認してみましょう。
以下は、結果を表示するためのプログラム例です。

c idct2.c

#include <stdio.h>
#include <math.h>
int main(void) {
    int i, j; // loop
    double coefficients[10] = {7.778175, -4.036036, -0.000000, -0.432302, -0.000000, -0.141421, -0.000000, -0.057185, 0.000000, -0.016036}; // DCT係数
    double resample[10]; // IDCTから生成される信号
    const int size = 10; // 配列サイズ
    
    // idct
    { /* 省略 */ }
    
    // 計算結果の表示
    printf("coefficients -> resample\n");
    for (i = 0; i < size; i++) {
        printf("%6f,\t%6f\n", coefficients[i], resample[i]);
    }
    
    return 0;
}

結果はこんな感じで表示されると思います。

coefficients -> resample
7.778175,	1.000001
-4.036036,	2.000000
-0.000000,	3.000000
-0.432302,	4.000001
-0.000000,	5.000000
-0.141421,	6.000001
-0.000000,	7.000000
-0.057185,	8.000000
0.000000,	9.000001
-0.016036,	10.000000
Program ended with exit code: 0

DCT, IDCTの関数化

実際にDCT, IDCTを使って電子透かし等を埋め込む際には、関数化したほうが、バグが発生しにく、プログラムの可読性がかなりあがります。
「関数がちょっと苦手だな」って思う場合でも、敬遠せずにここで練習すると思って関数を作ってみましょう。

DCTの関数化

DCT関数のプログラムは以下のようになります。

c dct.c

void dct(double signal[], double write_to[], const int size) {
    int k, m;
    
    if (signal == NULL || write_to == NULL || size <= 0) {
        fprintf(stderr, "line: %d, %s() error: The function argument contains an incorrect value.\n", __LINE__, __FUNCTION__);
        return;
    }
    
    write_to[0] = 0;
    for (m = 0; m < size; m++) {
        write_to[0] += signal[m];
    }
    write_to[0] *= sqrt(2) / size;
    
    for (k = 1; k < size; k++) {
        write_to[k] = 0;
        for (m = 0; m < size; m++) {
            double vector = (2 * m + 1) * k * M_PI / (2 * size);
            write_to[k] += signal[m] * cos(vector);
        }
        write_to[k] *= 2.0 / size;
    }
    return;
}

関数の解説を書いていきます。

引数, 返り値を設定し、関数を宣言する
void dct(double signal[], double write_to[], const int size)
関数を作る場合は、引数, 返り値を設定する必要があります。
引数は計算するのに何の値が必要なのか考えるのが重要です。（実際にプログラムを書きながら設計する場合は、必要なパラメータ（引数）があったら少しずつ書き加えていく感じで作っていきます。）

DCTの場合は、まず時間領域の信号を受け取る必要があります。
関数化する前は、x[]という名前にしてましたが、関数ではより具体的なsignal[]に名前を変更しました。

write_to[]は計算結果を書き込むための配列を受け取っています。
c言語の場合、関数は基本的に1つの値しか返すことができないため、計算結果が配列になる場合は、引数であらかじめ書き込む先の配列を受け取る必要があります。

const int sizeは、配列の大きさを受け取っています。
constをつけなくても良いのですが、一応、配列のサイズは不変であるため定数表記にしておきます。

頭についているvoidですが、これは、返り値の型を表しています。
通常であれば計算結果を返すことが多いのですが、今回の場合はwrite_to[]で計算結果を返してしまうので、他に返す値はありません。
返す値がない場合は、voidにすれば値がないというふうに関数を宣言することができます。
イレギュラな引数を排除する
```
if (signal == NULL || write_to == NULL || size <= 0) {
    fprintf(stderr, "line: %d, %s() error: The function argument contains an incorrect value.\n", __LINE__, __FUNCTION__);
    return;
}
```
ここでは、おかしな引数が入っていないか確認する作業を行っています。
正しい引数が入っていない場合は、計算ができないので何がおかしいのかを表示し、計算を中止する必要があります。

基本的にはないのですが、配列であるsignal, write_toは、うまくメモリが確保できない場合は、NULL値が入る場合があります。
メモリがうまく確保できていない場合も、結果を書き込むことはもちろん、データを読み出すこともできないのでif文でチェックします。

また、sizeが0以下であった場合は、配列のサイズが0以下と、あり得ない状態であるため、計算をすることができませんのでif文でチェックします。

最後にエラーメッセージを追加することで、どこで問題が起きたかを可視化することができるのでメッセージを追記しておきます。
1.で設定した引数を使ったプログラムに書き換える
1.で名称が変わった変数名がありますので、変更した名称にした名前にプログラムを変更し、関数内にプログラムを移しましょう。
上記に記載したdct()関数では、[1]に乗っ取りfor文で使う変数名も変更しています。
もし、自力でプログラムを書き換えている場合は注意してください。

IDCTの関数化

IDCT関数のプログラムは以下のようになります。

c idct.c

void idct(double coefficients[], double write_to[], const int size) {
    int k, m;
    
    if (coefficients == NULL || write_to == NULL || size <= 0) {
        fprintf(stderr, "line: %d, %s() error: The function argument contains an incorrect value.\n", __LINE__, __FUNCTION__);
        return;
    }
    
    double dc_value = coefficients[0] / sqrt(2);
    
    for (m = 0; m < size; m++) {
        write_to[m] = dc_value;
        for (k = 1; k < size; k++) {
            double vector = (2 * m + 1) * k * M_PI / (2 * size);
            write_to[m] += coefficients[k] * cos(vector);
        }
    }
    return;
}

意識して書いたところに関しては、dctとほとんど一緒なのでidct1.cのプログラムと違うところだけ説明します。

直流成分をあらかじめ計算しておく
dc_valueという変数をおき、あらかじめ計算しておくようにしました。
直流成分は、for文の中で毎回使われるのものではありますが、変数によって変化するものではありません。
for文の中で毎回計算してしまうとほんの少しではありますが、計算遅延が発生するので、変数をおき、あらかじめ計算して呼び出すような形にしました。

関数の利用

以下のコードは関数利用例です。

c dct_idct.c

#include <stdio.h>
#include <math.h>

void dct(double signal[], double write_to[], const int size) { /* 省略 */ }

void idct(double coefficients[], double write_to[], const int size) { /* 省略 */ }

int main(void) {
    double x[10] = {10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 15.0, 14.0, 13.0};
    double coefficients[10];
    double resample[10];
    int size = 10;
    
    // DCT
    dct(x, coefficients, size);
    
    // 透かしの埋め込み
    { /* 何らかの処理で周波数領域に透かし情報を埋め込む */ } // 今回は何にもしない
    
    // IDCT
    idct(coefficients, resample, size);
    
    printf("signal -> coefficients -> resample\n");
    for (int i = 0; i < size; i++) {
        printf("%6f,\t%6f,\t%6f\n", x[i], coefficients[i], resample[i]);
    }
    
    return 0;
}

上記のコードを見てみてどうでしょうか？
main関数でfor文を利用しないでDCT, IDCTができています。
短いコードで処理の流れがかなり見やすくなっていると思います。

電子透かしでは実際に透かしを埋め込む場合、音声データなどをブロック分割し、繰り返しDCTを行ったりすることもあると思います。
あらかじめ関数化しておけば、２重ループなどを意識せずにコードも書けますし、誤って使わない変数にアクセスするリスクも避けることができます。

こういった決まったものでなくても、ある程度決まった処理を繰り返し使う場合は、今回のように関数化することも是非考えてみてください。

まとめ

今回、離散コサイン変換 DCTについて調べ、プログラムを１から書いてまとめてみました。
詳しく調べてみると、DCTは１〜４系と、いくつか用途によって式が分かれているみたいです。
それぞれの式の特徴までは調べることはできませんでしたが、余裕があったら調べてみるのも良いですね。
ぜひ、この記事を読んで満足せず、英語ではありますが参考文献[1]も目を通してみてください。

この記事を書く上でところどころ、自己解釈が入っており、誤っているところがあるかもしれません。
もし見つけた場合は、遠慮せずコメントにてお知らせいただけると幸いです。

おまけ

TeXで使用できる式置き場

DCT

$X(0) =\frac{\sqrt{2}}{N} \sum\limits_{n=1}^{N-1} x(n) $
$X(k) =\frac{2}{N} \sum\limits_{n=1}^{N-1} x(n)\cos{\frac{(2n+1)k\pi}{2N}}, k=1,2,...,(N-1) $

X(0) =\frac{\sqrt{2}}{N} \sum\limits_{n=1}^{N-1} x(n)
X(k) =\frac{2}{N} \sum\limits_{n=1}^{N-1} x(n)\cos{\frac{(2n+1)k\pi}{2N}}, k=1,2,...,(N-1)

IDCT

$x(n) =\frac{1}{\sqrt{2}}X(0)+ \sum\limits_{k=1}^{N-1} X(k)\cos{\frac{(2n+1)k\pi}{2N}}, n=0,1,...,(N-1) $

x(n) =\frac{1}{\sqrt{2}}X(0)+ \sum\limits_{k=1}^{N-1} X(k)\cos{\frac{(2n+1)k\pi}{2N}}, n=0,1,...,(N-1)

参考文献

N. Ahmed, T. Natarajan and K. R. Rao, "Discrete Cosine Transform," in IEEE Transactions on Computers, vol. C-23, no. 1, pp. 90-93, Jan. 1974, doi: 10.1109/T-C.1974.223784.
Strang, Gilbert. "The discrete cosine transform." _SIAM review_41.1 (1999): 135-147.
安達丈晴, 長谷川まどか, 加藤茂夫, DCTを利用した静止画像の電子透かし法についての検討, 映像情報メディア学会技術報告, 1999, 23.62 巻, セッションID MIP99-70, p. 17-22, 公開日 2017/06/23, Online ISSN 2424-1970, Print ISSN 1342-6893, https://doi.org/10.11485/itetr.23.62.0\_17, https://www.jstage.jst.go.jp/article/itetr/23.62/0/23.62\_17/\_article/-char/ja