【離散コサイン変換】Part1. DCTをゆっくりじっくり理解する

はじめに

離散コサイン変換(以下DCTと呼びます)といえば画像処理や音声分野での応用が多く、計測や制御の現場で直接的に登場することはあまり無いように思います。

そのため、これまで「要するにフーリエ変換じゃないの？」とか、「よく分からないけどcosだけで信号を表現できるらしい！」程度の理解ですませていました。

ですが、圧縮センシングやスパース性を活かした信号処理を計測や制御に応用している事例を見聞きしていく内に、「思ってたより応用範囲広そう！」とか、「知らないとちょっと恥ずかしい、、、」と思ったので今回勉強しました。

基本的に自身の備忘のために書いた記事ですが、私と同じように「DCTが気になるけどよく知らない人」の参考にもなれば嬉しいです。

注）今回扱うDCTは全てタイプ2です。

離散フーリエ変換と離散コサイン変換の定義式

DCTも周波数領域で信号を分解するということなので、離散フーリエ変換(DFT)との違いを式から見ていきましょう。

とりあえず離散フーリエ変換がどんな式だったかな？というところから始めます。
離散フーリエ変換を含めたフーリエ変換ファミリーの定義式を並べてみます。
※以下、$f=\frac{1}{T}$, $ω=2πf$とします。

＜フーリエ級数展開＞　やってること：周期信号を三角関数で展開

・三角関数を用いた形式

$$f(t)=a_0+\sum_{n=1}^{∞}(a_ncos(nωt)+b_nsin(nωt))$$

・複素指数関数を用いた形式

$$f(t)=\sum_{n=-∞}^{∞}c_ne^{jnωt}$$

＜フーリエ変換＞　やってること：非周期信号を無限大の周期を持つ信号と見なしてフーリエ級数展開

・三角関数を用いた形式

$$X(ω)=\int_{-∞}^{∞}x(t)(cos(ωt)-jsin(ωt))dt$$

・複素指数関数を用いた形式

$$X(ω)=\int_{-∞}^{∞}x(t)e^{-jωt}dt$$

＜離散時間フーリエ変換(DTFT)＞　やってること：時間を離散的に扱ったフーリエ変換、周波数は連続

・三角関数を用いた形式

$$X(ω)=\sum_{n=-∞}^{∞}x[n](cos(nω)-jsin(nω))$$

・複素指数関数を用いた形式

$$X(ω)=\sum_{n=-∞}^{∞}x[n]e^{-jωn}$$

＜離散フーリエ変換(DFT)＞　やってること：時間に加えて周波数も離散的に扱ったフーリエ変換
　※信号は周期N
・三角関数を用いた形式

$$X[k]=\sum_{n=0}^{N-1}x[n](cos(\frac{2π}{N}kn)-jsin(\frac{2π}{N}kn))$$

　　　　　ωを使って表すと、

$$X[ω_k]=\sum_{n=0}^{N-1}x[n](cos(ω_kn)-jsin(ω_kn)) , ω_k=\frac{2π}{N}k (k=0,1,2,...)$$

・複素指数関数を用いた形式

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j\frac{2π}{N}}{kn}$$

　　　　　ωを使って表すと、

$$X[ω_k]=\sum_{n=0}^{N-1}x[n]e^{-jω_kn} , ω_k=\frac{2π}{N}k (k=0,1,2,...)$$

離散フーリエ変換含め、フーリエ変換はとにかくcos/sin両方を使った式になりますね。
ではDCTはどうでしょうか？DFTと比較してみましょう。
※周期Nの信号に対しての式です。

＜離散コサイン変換(DCT) タイプ2＞

$$X[k]=\sum_{n=0}^{N-1}x[n]cos(\frac{π}{N}k(n+\frac{1}{2}))$$

　※DFT再掲
　　

$$X[k]=\sum_{n=0}^{N-1}x[n](cos(\frac{2π}{N}kn)-jsin(\frac{2π}{N}kn))$$

コサイン変換なので当たり前ですが見ての通り、cosしか出てきませんね。
ということは、周波数領域で信号を分解することは変わらないが、複素数ではなく実数のみで信号を表せるように分解しているということになりますね。

DCTの定義式導出

DFTの式を元に上記DCT(タイプ2)の定義式を導いてみます。

・STEP1：長さNの信号$x[n]$を拡張して、偶関数である長さ2Nの信号$\tilde{x}[n]$を作ります。
　　　　こんな感じですね。

$$\tilde{x}[n] =\begin{cases}x[n] & \text{for } 0 ≦ n < N \\x[2N-1-n] & \text{for } N ≦ n < 2N\end{cases}$$

・STEP2：拡張した信号に対して2N点のDFTを行います。

$$X_k=\sum_{n=0}^{2N-1}\tilde{x}[n]\cdot e^{-j\frac{2π}{2N}kn}, k=0,1,...,2N-1$$

ここでポイントですが、2Nまで拡張して偶関数になった状態の信号に対して2N点のDFTを行うことで、結果的に虚部が打ち消し合って実部だけが残ります。イメージ的には下記図のようになります。

・STEP3：DFT結果を三角関数を用いた形式に書き換えます。

$$X_k=\sum_{n=0}^{2N-1}\tilde{x}[n]\cdot (cos(\frac{2π}{2N}kn) - jsin(\frac{2π}{2N}kn)), k=0,1,...,2N-1$$

虚部は打ち消し合うのでsinの項は消えて、ついでにcos内を約分します。

$$X_k=\sum_{n=0}^{2N-1}\tilde{x}[n]\cdot (cos(\frac{π}{N}kn)), k=0,1,...,2N-1$$

ここまでで、元信号を拡張して2N点の偶関数にした信号のDFTをcosのみで表している状態ですね。

・STEP4：次の目標はSTEP3までの式をN点の長さで、$x[n]$を使った形で表すことです。

まずSTEP3の式を前半と後半に分けて、次のように書き直します。

$$X_k=\sum_{n=0}^{N-1}\tilde{x}[n]\cdot (cos(\frac{π}{N}kn))+\sum_{n=N}^{2N-1}\tilde{x}[n]\cdot (cos(\frac{π}{N}kn))$$

$\tilde{x}[n]$は偶関数なので、$\tilde{x}[n]=x[n] (0≦n<N)$と$\tilde{x}[n]=x[2N-1-n] (N≦n<2N)$が成り立ちます。
なので、$m=2N-1-n$と置くと、$n=2N-1-m (0≦m≦N-1)$を使って、後半の項は

$$\sum_{n=N}^{2N-1}\tilde{x}[n]\cdot (cos(\frac{π}{N}kn))=\sum_{m=0}^{N-1}x[m]\cdot (cos(\frac{π}{N}k(2N-1-m)))$$

と書けます。

そして、cosは$2π$周期のため

$$cos(\frac{π}{N}k(2N-1-m))=cos(2πk-\frac{π}{N}k(m+1))=cos(\frac{π}{N}k(m+1))$$

と表せて、前半部分の項($\tilde{x}[n]をx[n]$に書き直しています)と合わせると、

$$X_k=\sum_{n=0}^{N-1}x[n]\cdot (cos(\frac{π}{N}kn))+\sum_{m=0}^{N-1}x[m]\cdot (cos(\frac{π}{N}k(m+1))$$

となります。

$0≦n≦N-1$と$0≦m≦N-1$の範囲が同じなので、変数を$n$に揃えて

$$X_k=\sum_{n=0}^{N-1}x[n]\cdot (cos(\frac{π}{N}kn)+cos(\frac{π}{N}k(n+1)))$$

そしてcosの和の公式から$\sum x[n]$以降は、

$$cos(\frac{π}{N}kn)+cos(\frac{π}{N}k(n+1))=2cos(\frac{π}{N}k(n+\frac{1}{2}))\cdot cos(\frac{π}{2N}k)$$

と書き直せるので、$n$に依存しない定数$2$と$cos(\frac{π}{2N}k)$を外に出すと、

$$X_k=2cos(\frac{π}{2N}k)\cdot \sum_{n=0}^{N-1}x[n]\cdot cos(\frac{π}{N}k(n+\frac{1}{2}))$$

となって、定数部分$2cos(\frac{π}{2N}k)$を除いたところがDCTタイプ2の定義式となります。長かった、、、

数式から見たDCTの特徴

一般的に言われるDCTの性質というかメリットは以下のようなものだと思います。
・エネルギーが低周波に集中しやすい
・情報の圧縮に適している

これらのよく言われる特徴を下記DCTの定義式から導き出してみましょう。
再掲：DCTタイプ2の定義式

$$X_k=\sum_{n=0}^{N-1}x[n]cos(\frac{π}{N}k(n+\frac{1}{2}))$$

＜右辺の意味＞

まず右辺を考えてみます。右辺の意味を日本語で書き下すと下記のようになりますよね。
「元の信号$x[n]$と$cos(\frac{π}{N}k(n+\frac{1}{2}))$の内積を取る」

これを言い換えると、下記のようになります。
「元の信号$x[n]$と$cos(\frac{π}{N}k(n+\frac{1}{2}))$の類似度を測っている」

この時、$cos(\frac{π}{N}k(n+\frac{1}{2}))$の方はkの値が大きくなるほど周期の短い、高周波なcos波になっていきますよね。

一方で元の信号$x[n]$は変わらないので、$x[n]$が直線だったり周期の長い正弦波の場合はkの値が小さいcos波との内積が大きくなり、逆に$x[n]$が周期の短い、高周波成分を多く含む信号だった場合はkの値が大きいcos波との内積が大きくなります。

それを図にすると下記のようになります。
k=0とk=1の図ですが、上段(低周波)から中段、下段(高周波)になるほど元信号の周波数は高くなっていますが、同時に内積は小さくなっていますね。

これがDFTの場合はcosと同時にsinとの内積も取ることになるので、元信号との類似度はもう少し複雑になりますが、DCTの場合はcosだけなので単純に考えることができますね。

＜低周波へのエネルギー集中＞

ここでDCTタイプ2の定義式を導出した時の操作を思い出すと、最初に長さNの信号を拡張して長さ2Nの偶関数になった信号を作りましたよね。

DCTの場合、この偶関数への拡張操作によって信号の端と端がきれいに繋がって滑らかになり、結果として対象の信号は高周波成分が抑えられた構造になります。

そこで先程の話に戻りますが、元の信号が高周波成分を多く含むほどkの値が大きいcos波との内積が大きくなり、低周波成分がほとんどの場合はその逆のことが起こるのでした。

ということは、高周波成分が抑えられた構造の信号においては、kの値が小さいcos波との内積が大きくなりやすく、kの値が大きいcos波との内積は小さいかほぼゼロになりやすいということです。

そして内積を求めた結果$X_k$はDCT係数ですから、$X_k$の絶対値を二乗して全てのkについて総和を取ると、信号全体のエネルギーになりますよね(パーセバルの等式ですね)。

しかし、そもそもkの値が大きいcos波との内積はごく小さい、またはゼロなので、そのときの$X_k$も小さくなり、よってそれの二乗によって求められる、kの値が大きいすなわち高周波成分のエネルギーが信号全体のエネルギーに占める割合も小さくなります。

ここまでがDCTの特徴：「エネルギーが低周波に集中しやすい」ことの説明になります。

＜DCT係数のスパース性＞

低周波へのエネルギー集中と地続きの話になりますが、滑らかな信号の場合、kの値が大きいcos波との内積が(ほぼ)ゼロになる、という説明をしました。

ということはつまり、繰り返しになりますが元信号に占めるその(高い)周波数のエネルギーの割合がごく小さいということで、そうするとその分は元信号の情報として省いてもさほど問題無いわけです。

すなわち、実質的に同じ信号をより少ない非ゼロのDCT係数でほぼ正確に表現できるということで、
ひとことでまとめると、「滑らかな信号に関するDCT係数はスパースになる」ということです。

同じものをより少ない情報量で表現できるというのは、とりもなおさず情報の圧縮そのものですよね。

これがDCTの特徴：「情報の圧縮に適している」ことの説明になります。

最後に、元信号$x[n]$の違いとkの値による内積の変化をk=0~7まで見てみましょう。
上段(低周波)の場合、kの値が大きくなるにつれて内積が顕著に小さくなっていくのが分かりますね。

まとめ

ここまで、DFTを元にして離散コサイン変換(タイプ2)の定義式を導出し、数式から見たDCTの性質を考えてきました。

個人的な感想ですが、巷で言われているDCTの特徴について割と腑に落ちたような感じで、ゆっくりじっくり、ある程度理解できたと思います。

ただ、信号処理はやっぱり応用してナンボなところがあるので、次回記事では結局DCTの特徴が分かったからって何に使えて、何が嬉しいの？というところを見ていきたいと思います。