DCT(離散コサイン変換) - abwiki @ ウィキ (ActiveBasic非公式wiki)

離散コサイン変換(DCT)とは、波はSIN関数の合成で表現できるというフーリエ級数の
考え方に従った変換で、物理学的な意味は時間変動する波から周波数成分へと変換するための方法。

同じ手法として比較的良く知られているのは離散フーリエ変換(DFT)や高速フーリエ変換(FFT)
であるが、離散コサイン変換はコンパクトで実数のみで計算ができることから実装しやすく、
JPEGやMPEGなどの動画、MP3などの音声圧縮などに多く用いられている。
MMXなどで計算することも可能で、コーデックのDCTがMMX/SSEに対応している場合は、
PentiumMMXなどのCPUでもMPEG1動画を容易に再生することができた。

離散コサイン変換(DCT)は時間変動する波形を周波数成分に変換し、逆離散コサイン変換
(InverseDCT)は、周波数成分の分布から時間波形を再構成する。
時間波形を周波数変換する場合は、フーリエ級数の考え方に従い、さまざまな周波数の
振動波形(DCTの場合はCos関数)にたいして源信号との間で積分し、その後に周波数成分を
得ている。
通常、DCTと言う場合はDCT-2の式を意味し、IDCTがDCT-3の式を表す。下記のDCTはDCT-2となり、
DCTはデータ数は任意に指定できる。DCTで得られた数列に大して絶対値/Abs/を計算する
ことでパワースペクトルを得られる。
DCTは複素関数を含まないので、コードは比較的短く簡単。

高速DCTは昔からいくつかの方法が提案されテいるが、
高速化の方法としては、FFTのようにバタフライ演算/バタフライソート演算/を用いる方法、
単純に求めたいデータを間引いて計算を少なくする方法、MPEGで用いられているような
整数DCTを用いる方法、スレッドやparallel~forなどを用いる方法などが考えられる。

DCT

#console

'DCT
'離散コサイン変換サンプル

Dim dat(10) As Double
Dim ret(10) As Double
Dim a As Double
Dim b As Double
Dim d As Double
Dim pi As Double

Dim i As Long
Dim j As Long
Dim L As Long

'定数
pi=3.14159
L=10

'データ定義(時間波形)
dat(0)=1
dat(1)=1
dat(2)=1
dat(3)=0
dat(4)=0
dat(5)=0
dat(6)=1
dat(7)=1
dat(8)=1
dat(9)=0


'--- DCT Main
a=2*L
d=Sqr(2/L)

for i=0 to L
b=i*pi
ret(i)=0

for j=0 to L
ret(i)=ret(i)+dat(j)*Cos( ((2*j+1)*b)/a )
Next j
ret(i)=ret(i)*d

if i=0 then ret(i)=ret(i)/Sqr(2)
Next i

'結果表示(周波数成分)
for i=0 to L
print ret(i),
Next i

iDCT

#console

'iDCT
'逆離散コサイン変換サンプル
'

Dim dat(10) As Double
Dim ret(10) As Double
Dim a As Double
Dim b As Double
Dim d As Double
Dim pi As Double

Dim i As Long
Dim j As Long
Dim L As Long

'定数
pi=3.14159 'PI
L=10 'Array Length(データ数)

'入力データ定義(周波数成分)
dat(1)=1

'--- idct main
a=2*L
d=Sqr(2/L)

for i=0 to L
b=2*i+1
ret(i)=0

for j=0 to L
ret(i)=ret(i)+dat(j)*Cos( (j*pi*b)/a )
if i=0 then ret(i)=ret(i)/Sqr(2)
Next j
ret(i)=ret(i)*d
Next i

'結果表示(時間波形)
for i=0 to L
print ret(i),
Next i

• バタフライソート演算のAppropriateなせつめい。

DCTやDFTを高速化する際に行なわれる計算で、計算式のうち複数の同じような値の式を
省くために行なわれる。

時間波形を示すデータに対してDCT/DFTを行なう場合、DCTでは元のデータに
対してcos(wt)の項で積を求めるが、この際に幾つか異なる計算式のなかに同じ項
(factor)が存在していることがある。

下の適当な例では(式は適当なので信用してはいけない。あくまでもイメージをつかめればよい)、
データ数n=4のDCTの場合、DCT-2の式のΣを丁寧に展開する訳だが、BASICプログラムの場合は
for文の書かれているループ項だ。

for j=0 to L
ret(i)=ret(i)+dat(j)*Cos( ((2*j+1)*b)/a )
Next j

ここを詳しく見るためにforを使わずに書けば、下記のような計算と等しい.

F(0)=dat(0)*cos(w0t0) + dat(1)*cos(w0t1) + dat(2)*cos(w0t2) + dat(3)*cos(w0t3)
F(1)=dat(0)*cos(w0t0) + dat(1)*cos(w1t1) + dat(2)*cos(w2t2) + dat(3)*cos(w3t3)
F(2)=dat(0)*cos(w0t0) + dat(1)*cos(w2t1) + dat(2)*cos(w4t2) + dat(3)*cos(w6t3)
F(3)=dat(0)*cos(w0t0) + dat(1)*cos(w3t1) + dat(2)*cos(w5t2) + dat(3)*cos(w7t3)

縦のF(n)は求められるスペクトル数列であり、横のdat(n)は時間波形の
サンプルデータ、そしてデータに対してCos関数等が書かれているツイストファクターが存在する。
このなかで幾つかのFactorに同じような計算が現れることがあり、それらの計算を省く
事が出来れば、必要な計算処理数が小さくなる。sin/cosは周期関数だから共通の項があっても
おかしくは無い。

実際に同じ項が存在するので、それを省いてF(n)を求める事が出来れば計算量を少なく
することができる。

この特性を詳しく見てみると、Cos(wt)等の同じ項の計算は、別の周波数F(n)を求める式の項の中に含まれていたりする、
つまりF(0)を求める際にCos(wt)等のツイストファクタの項はF(3)などに現れていたりするという事。

なので同じ計算はたすきがけのようなクロスした線で結ばれた位置にあり、それらを参照すれば
重複する計算を削減でき計算量が減る。

しかしF(0)とF(3)では計算順序が前後関係として前と後ろにあって、求めたい値がまだ計算
されていないことになる。これは困ったということで、それなら共通の項が現れる法則に
したがって、初めからソートしておいて、その後に計算すればいいという発想が生まれる。

このクロスした関係を図に表すとバタフライのように線が書かれる事からバタフライ演算と呼ぶ。
重複した計算を省くバタフライ演算の際は、実際にはソートなどを行なって計算式(計算順)を変える。
これをビットソートと呼び、偶数と奇数に分けビット列を参照しながらソートする。
ビットソートはバタフライ演算で計算を短くする際に用いられる前処理となる。

バタフライ演算の説明で良く示されるクロス線の書かれた図は、上の代数式の各項の関係を線で
結んだもので、どの項とどの項が関係しているかを線のクロスで図示してある。
(実際の図は、上の式を多少変形しているのでクロス線に書かれるのだが、上の式のままでは多分
変な線の図になると思う)

本来は行列式で表したほうがすっきりとするが、バタフライという名前の由来を知るため
には行列式を使わない表現のほうが意味が判りやすいだろう。

このバタフライソートによって計算が省かれるので、例えばn=4の場合

F(n)=dat(0)*cos(w0t0) + dat(1)*cos(w2t1) + dat(2)*cos(w4t2)

等というように本来は4回の所を3回程度に計算式を短縮できる。(この式も適当だが)
n=8の場合も同様に3-4回に短縮される。もしバタフライソート演算をしなければ8回計算しなくてはならない事になる。単純に8*8(8^2)という計算量となる
だろう事が予想される。一般にDCT/DFTがn^2という計算量を持つという理由だ。

バタフライソート演算では重複した計算を省くので8*3相当の式程度に計算量は短くなる。大体のイメージはわかっただろうか？
書いている本人も十分に分かっていないのだが、多分合っている。