高速ハートレー変換に基づく高速なDCTアルゴリズム

第26回信号処理シンポジウム 2011年11月16日∼18日(札幌)
高速ハートレー変換に基づく高速な DCT アルゴリズム
Fast Hartley Transform based Fast DCT Algorithm
鈴木大三 1
田中雄一 2
池原雅章 3
阿曽弘具 1
1：日本大学工学部電気電子工学科
2：宇都宮大学工学部情報工学科
3：慶應義塾大学理工学部電子工学科
Taizo SUZUKI1
Yuichi TANAKA2
Masaaki IKEHARA3
Hirotomo ASO1
1: Department of Electrical and Electronics Engineering, College of Engneering, Nihon University
2: Department of Information Science, Utsunomiya University
3: Department of Electronics and Electrical Engineering, Keio University
アブストラクト
離散コサイン変換（DCT: discrete cosine transform）は効率の良い周波数変換として知られて
案する．そして，回転行列の特性を用いて，その構造の冗
おり，ソフトウェア／ハードウェアで高速処理には，乗
長性を除去した後，乗算器を使用しないリフティング構
算器を使用しない構造が望ましい．本論文は，ロッシー
造 [8] に近似することで，高速な DCT を得る．その DCT
（非可逆）画像符号化のための乗算器を使用しない高速な
は，乗算器を使用せず加算器とビットシフト器（以下，シ
DCT の実現を提案する．まず高速ハートレー変換（FHT:
fast Hartley transform）に基づいた新たな DCT を提案
フト器）のみで構成することができる．提案 DCT は，高
する．次に回転行列の特性を利用して，その構造の冗長
想的なオリジナル DCT の持つ特長を継承し，IntDCT [5]
性を除去する．そして回転行列を乗算器を使用しないリ
や BinDCT [6] のような従来の高速アルゴリズムよりも
フティング構造に近似することで，乗算器を使用しない
良好な符号化性能を示す．
DCT が実現される．最後に従来の DCT との画像符号化
比較によって，提案 DCT の有用性を実証する．
1
transform）[9] に注目し，それに基づく新たな DCT を提
い符号化利得，DC 漏れなし，対称な基底関数といった理
表記： I，J，·T はそれぞれ，単位行列，反転行列，行
列の転置を示す．
はじめに
離散コサイン変換（DCT：discrete cosine transform）[1]
はエネルギー圧縮に関して十分な性能を示すことができ，
基礎理論
2
2.1
またソフトウェア／ハードウェア実装のために多くの高
離散コサイン変換（DCT：discrete cosine
transform）
速なアルゴリズムが提案されている．その結果，DCT は
画像／動画像処理やその他の様々な分野で応用されてい
DCT [1] にはいくつかのタイプがあるが，本論文ではタ
る．JPEG，MPEG，H.26x シリーズ [2–4] にも採用され
イプ II（DCT-II）とタイプ III（DCT-III）についてのみ
ており，そのような多くの国際標準においても中核を成
記述する．これらはそれぞれ一般的に DCT 及び逆 DCT
す技術となっている．
と呼ばれる変換であり，JPEG，MPEG，H.26x シリー
通常の従来の DCT は，整数の時間信号を実数の周波数
ズ [2–4] のような画像／動画像符号化標準にも採用されて
信号へ置き換える変換である．そのため，ソフトウェア
いる．M 分割 DCT 行列及び逆 DCT 行列をそれぞれ C
／ハードウェアで実装される際の計算コストや消費電力， 及び D とするならば，その m 行 n 列要素は以下のよう
特に実数で表される乗算器を軽視することはできない．特
に表される．
に携帯機器を実装する際などに，消費電力の問題は重要
な項目となる．つまり，実数で表される乗算器を可能な
限り削減，そして完全に除去できるような DCT の新しい
アルゴリズム [5–8] を創出することは重要な課題である．
本論文では，高速ハートレー変換（FHT：fast Hartley
- 584 -
√
[C]m,n
[D]m,n
(
)
2
m (n + 1/2) π
cm cos
=
M
M
√
(
)
2
(m + 1/2) nπ
cn cos
=
M
M
図 1: リフティング構造（白丸：ラウンディング処理）
図 2: 2 進係数の乗算器とラウンディング処理を持つリフ
ティング構造の加算器とシフト器への近似（≫ n: n ビッ
トシフト器）
−1
ただし，D = C = C ，0 ≤ m, n ≤ M − 1，


 √1 (n = 0)
 √1 (m = 0)
2
2
及び cn =
cm =

1
(m ̸= 0)
1
(n ̸= 0)
T
3 ビットシフト器，4 ビットシフト器，7 ビットシフト器
の加算に置換できる．図 2 にその様子を示す．このよう
である．単純のために，M = 2n (n ∈ N) と定義する．ま
な近似が行われても，リフティング構造自体の完全再構
成は保持されていることに注意されたい．
た，DCT 及び逆 DCT の高速アルゴリズムが多数報告さ
れている [10–12]．
2.2
高速ハートレー変換（FHT: fast Hartley trans-
3
form）に基づく乗算器を使用しない高速な DCT
乗算器を使用しないリフティング構造
JPEG2000 [13] や JPEG XR [14] を代表とする次世代
本章では，ロッシー（非可逆）画像／動画像符号化の
画像符号化標準における周波数変換では，しばしばリフ
ための乗算器を使用しない高速な DCT を，以下のような
ティング構造 [15] が用いられている．整数入力−整数出
手順で導出する．
力変換（整数変換）を実現し，ロスレス（可逆）モード
1. FHT に基づく新たな DCT の構造を示す．
を実現することのできる手法として，また，高速処理を
実現する構造として有用なためである．図 1 にリフティ
2. 回転行列の特性を利用して，その DCT を構成する
冗長な回転行列を除去する．
ング構造の基本構造の分割側及び合成側の概略図を示す．
この場合，以下のような式で表すことができる．
yj (n) = xj (n), yi (n) = xi (n) + round{P xj (n)}
3. いくつかの回転行列を乗算器を使用しないリフティ
ング構造に近似し，全てのスケーリング係数を量子
化部に移動することで，乗算器を使用しない高速な
zj (n) = yj (n), zi (n) = yi (n) − round{P yj (n)}
このとき，P 及び round{·} はそれぞれリフティング係数
DCT を実現する．
及びラウンディング処理（丸め処理）を表す．また，この
場合のリフティング行列及びその逆行列は，それぞれ以
3.1
FHT に基づく DCT
下のように表すことができる．
[
]
1 P
0 1
[
1
0
及び
P
1
]−1
FHT [9] はディジタル信号処理のアルゴリズムやシス
[
]
1 −P
=
0
1
テムの解析，設計，処理などに用いられる変換の一つであ
る．対称直交行列で表される M 分割 FHT 行列 H（H−1 =
実数のリフティング係数の乗算器を k/2n （k, n ∈ N）
HT = H）の m 行 n 列要素は以下のように表される．
(
)
1
2mnπ
[H]m,n = √
cas
M
M
(
(
)
(
))
1
2mnπ
2mnπ
=√
cos
+ sin
M
M
M
のような 2 進係数に近似することが，高速な処理のため
の第一歩である．2 進係数 k/2n の乗算器とラウンディン
グ処理を持つリフティング構造は，加算器とシフト器に
近似される [8]．その結果，リアルタイムなソフトウェア
／エンコーダの実行の高速化や，回路規模の縮小化を図
ることができる．たとえば，リフティング係数 25/128 は，
4
3
ただし，このままでは出力信号が周波数帯域の順に並ん
でいないため，正しい周波数帯域の順に並べ替えた FHT
行列を Ĥ とし，以降これを使用する．また M = 8 の場
0
2 +2 +2
1
1
1
25
=
= 3 + 4 + 7.
128
27
2
2
2
合，FHT は図 3 に示されるように，回転角 π/4 を持つ 13
のように分解することができる．つまり，係数 25/128 の
乗算器とラウンディング処理を持つリフティング構造は，
個の回転行列によって単純化することができる．
今まで提案されていない DCT の新たな構造として，
- 585 -
図 3: 8 分割 FHT Ĥ.
↓
FHT に基づく構造を以下に示す．


1
0
0 0


0 JCJ 0 JS T
 Ĥ

C=

0
0
1
0


0
−SJ
0
(1)
↓
C
ただし，
(
[C]k,k = cos
(k + 1)π
2M
)
(
及び [S]k,k = sin
(k + 1)π
2M
)
であり，また 0 ≤ k ≤ M/2 − 2 である．つまり，DCT は
FHT 及び (M/2 − 1) 個の回転行列のみで構成することが
できる．M = 8 の場合の構造を図 4 上段に示す．
3.2
図 4: 提案する 8 分割 DCT：
（上段）FHT に基づく DCT，
（中段）冗長な回転行列の除去，
（下段）乗算器を使用しな
冗長な回転行列の除去
式 (1) の FHT に基づく DCT はいくつかの冗長な回転
い構造への近似
行列を含んでいることに注意されたい．M = 8 の場合，
以下のような回転行列の特性を利用して，図 4 中段に破
c(θ) = cos θ，s(θ) = sin θ，及び t(θ) = tan θ である．ま
線で示された 3 つの回転行列が除去される．
た，図 4 中段にある太線によって描かれた回転角 π/4 を
Θπ4 Θπ4 = I 及び
持つ回転行列にだけは，スケーリングの関係で，以下の
diag{1, −1}Θπ8 Θπ4 = Θπ8
ような一般的な回転行列のリフティング分解を適用する．
[
] [
][
][
]
1
0 1 s(θ)
1
0
c(θ) s(θ)
= 1−c(θ)
c(θ)−1
1 0 −1
1
s(θ) −c(θ)
s(θ)
s(θ)
ただし，
Θab
[
]
cos ( ab )
sin ( ab )
=
sin ( ab ) − cos ( ab )
この構造を図 5 下段に示す．そのようにして得られた実
である．このようにして冗長性が除去された構造は，回
転角 π/4，3π/16，π/8，π/16 を持つ 13 個の回転行列で
数のリフティング係数を k/2n のような 2 進係数に近似す
る．そして，図 4 中段の回転角の記載のない回転行列に
関しては，
構成される．
[
]
[ ( )
( )]
c π4
s π4
1 1 1
3.3 乗算器を使用しない構造への近似
( )
( ) =√
(3)
2 1 −1
s π4
−c π4
乗算器を使用しない DCT は以下の手順で得ることがで
√
きる．まず，図 4 中段において回転角の記載のある回転 のように，スケーリング係数 1/ 2 と加算及び減算に書
き換えることができる．最後に，式 (2) 及び (3) における
行列は，
√
][
][
[
] スケーリング係数 diag{c(θ), 1/c(θ)} 及び 1/ 2 を，それ
] [
c(θ)
0
1
0 1 t(θ)
c(θ) s(θ)
ぞれ量子化部に移動する．つまり提案 DCT は，図 4 下段
=
1
0
−s(θ)c(θ) 1 0
1
−s(θ) c(θ)
c(θ)
及び表 1 に示すように，乗算器に頼らず，加算器とシフ
(2)
のようにスケーリングと 2 ステップのリフティング構造
に分解される [6]．この構造を図 5 上段に示す．ただし，
ト器のみで実装することができる．1 また，表 1 における
1 ただし，[6] にしたがって本論文でも，実際は量子化ステップ及び
整数への丸め込みで取り扱われるべきスケーリング係数を，実数係数の
まま扱っている．
- 586 -
表 2: 8 分割 DCT における乗算器，加算器，シフト器の
数と符号化利得 Cg の比較
DCT
IntDCT
BinDCT
提案
[10]
[5]
[6]
DCT
加算器
52
26
0
45
0
40
0
38
シフト器
0
17
24
22
Cg
8.8259
8.7240
8.8244
8.8250
乗算器
図 5: 回転行列のリフティング分解：
（上段）スケーリン
うに定義される [16]．
グと 2 ステップのリフティング分解，
（下段）一般的なリ
Cg = 10 log10 (
∏M −1
フティング分解
k=0
σx2
σx2i ∥ fi ∥2
) M1
ただし，σx2 ，σx2i ，及び ∥ fi ∥2 はそれぞれ，入力信号の
分散，i 番目のサブバンド信号の分散，及び i 番目の合成
表 1: 提案 DCT のリフティング係数
実数係数
フィルタの基底関数のノルムを表している．表 2 に，8 分
2 進係数
p0
tan ( 3π
16 )
1
21
u0
3π
− sin ( 3π
16 ) cos ( 16 )
p1
tan ( π8 )
1
22
u1
− sin ( π8 ) cos ( π8 )
− 211
p2
π
)
tan ( 16
u2
π
π
− sin ( 16
) cos ( 16
)
p3
cos(π/4)−1
sin(π/4)
− 211 +
u3
sin ( π4 )
1−
p4
1−cos(π/4)
sin(π/4)
+
− 211 +
+
1
23
+
1
23
1
21
1
24
−
わらず，従来法である IntDCT [5] 及び BinDCT [6]（係
1
25
+
数は [7] におけるタイプ C1 を使用）よりも高い符号化利
1
25
得を示した．また，オリジナルの実数係数を持つ DCT に
は Chen の分解 [10] を用いた．2
1
24
− 212 +
1
22
利得 Cg の比較を示す．提案 DCT は少ない加算器にも関
1
25
+
1
23
+
割 DCT の場合の乗算器，加算器，シフト器の数と符号化
1
23
また図 6 に，オリジナル DCT 及び提案 DCT の周波数
1
24
+
−
特性の比較を示す．画像圧縮に重要な特性であるレギュ
1
25
ラリティ[16] を損なうことなく実現できていることが分
1
24
かる．
1
24
4.2
ロッシー画像符号化
本節では，PSNR によるロッシー画像符号化の比較に
おいて提案 DCT の優位性を示す．PSNR は以下の式で表
2 進係数は経験的に決定した．
される．
(
シミュレーション結果
4
4.1
PSNR [dB] = 10 log10
符号化利得と周波数特性
符号化利得は，圧縮アプリケーションに置いて考慮され
るべき最も重要な要因の一つである．高い符号化利得を持
つ変換は周波数分解能が高く，より多くのエネルギーを一
点に集中させることができ，結果として，ピーク信号雑音
比（PSNR: peak signal-to-noise ratio）のような客観的指
標で良好な性能を示す．DCT の符号化利得は最適な変換
であるカルーネン・レーヴェ変換（KLT: Karhunen-Loéve
transform）に近似されるため，乗算器を使用しない DCT
の符号化利得もオリジナルの DCT とできるだけ同等の符
号化利得を持つことが望ましい．符号化利得は以下のよ
2552
MSE
)
MSE は平均自乗誤差（Mean squared error）である．テ
スト画像には，Lena を代表とする 512 × 512 サイズの 8
ビット・グレースケール画像を用いた．また変換画像を符
号化するために，画像符号化のシミュレーションによく
用いられる SPIHT（The set partitioning in hierarchical
trees）[17] を採用した．
表 3 において，PSNR によるロッシー画像符号化の比
較を示す．また図 7 に，圧縮比 1:32 のときの Lena の拡
大画像の比較を示す．特に低ビットレートにおいて，提案
2 [8] はロッシー・ロスレス統合画像符号化を考慮した DCT である
ため，ロッシー画像符号化に特化した本論文の DCT との比較は行わな
かった．
- 587 -
図 6: 8 分割 DCT の周波数応答：
（左）DCT [10]，
（右）提案 DCT
表 3: ロッシー画像符号化比較（PSNR [dB]）
テスト画像
圧縮比
DCT [10]
IntDCT [5]
BinDCT [6]
提案 DCT
Barbara
1:64
1:32
1:16
24.45
26.95
30.70
24.31
26.68
29.96
24.44
26.88
30.45
24.42
26.88
30.47
Boat
1:64
1:32
1:16
26.06
28.66
31.91
25.90
28.43
31.51
26.02
28.58
31.55
26.03
28.60
31.64
1:64
1:32
27.36
29.38
27.16
29.17
27.31
29.25
27.39
29.31
1:16
31.99
31.66
31.60
31.69
1:64
1:32
28.71
31.91
28.46
31.48
28.65
31.69
28.69
31.78
1:16
35.64
34.83
34.86
35.04
1:64
1:32
28.22
31.43
28.02
31.16
28.18
31.27
28.19
31.31
1:16
34.49
34.00
33.92
34.04
Goldhill
Lena
Pepper
DCT は従来法と比べて視覚的にも同等以上の結果を示し
ていると言える．
5
まとめ
参考文献
[1] K.R. Rao and P. Yip, Discrete Cosine Transform
Algorithms, Academic Press, 1990.
本論文は離散コサイン変換（DCT: discrete cosine
[2] ISO/IEC 10918-1, Information technology - Digi-
transform）の乗算器を使用しないより高速なアルゴリズ
ムを提案した．はじめに，高速ハートレー変換（FHT: fast
tal compression and coding of continuous-tone still
images: Requirements and Guidelines.
Hartley transform）による新たな DCT の構造を紹介し
[3] ISO/IEC 14496-2, Information technology - Coding
た．次に回転行列の特性を用いて，その構造における全て
の冗長な回転行列を除去した．そして，いくつかの回転行
列を乗算器を使用しないリフティング構造へと置き換え，
またスケーリング係数を量子化部へ移動することにより，
of audio-visual objects - Part 2: Visual.
[4] ISO/IEC 14496-10, Information technology - Coding of audio-visual objects - Part 10: Advanced
Video Coding.
提案 DCT を得た．結果として，乗算器を用いず，またよ
り少ない加算器とシフト器によって構成しているにも関わ
らず，特に低ビットレートの画像符号化シミュレーション
[5] Y.J. Chen, S. Oraintara, and T.Q. Nguyen, “Video
compression using integer DCT,” Proc. of ICIP’00,
において，従来法と同等以上の結果を得ることができた．
- 588 -
Vancouver, British Columbia, Canada, Sep. 2000.
図 7: ロッシー画像符号化比較（テスト画像：Lena，圧縮比：1:32）
：
（上段）原画像，DCT [10]，IntDCT [5]，
（下段）
BinDCT [6]，提案 DCT
[6] J. Liang and T.D. Tran, “Fast multiplierless
[11] B.G. Lee, “A new algorithm to compute the discrete
approximations of the DCT with the lifting
scheme,” IEEE Trans. Signal Process., vol.49,
cosine transform,” IEEE Trans. Acoust., Speech,
Signal Process., vol.32, no.6, pp.1243–1245, June
no.12, pp.3032–3044, Dec. 2001.
1984.
[7] Y.J. Chen, S. Oraintara, T.D. Tran, K. Amaratunga, and T.Q. Nguyen, “Multiplierless approxi-
[12] Z. Wang, “On computing the discrete Fourier and
cosine transforms,” IEEE Trans. Acoust., Speech,
mation of transforms with adder constraint,” IEEE
Signal Process. Lett., vol.9, no.11, pp.344–347, Nov.
2002.
Signal Process., vol.33, no.4, pp.1341–1344, Apr.
1985.
[8] T. Suzuki and M. Ikehara, “Integer discrete cosine transform via lossless Walsh-Hadamard transform with structural regularity for low-bit-wordlength,” IEICE Trans. Fundamentals, vol.E93-A,
no.4, pp.734–741, Apr. 2010.
[13] ISO/IEC 15444-1, Information technology - JPEG
2000 image coding system: Core coding system.
[14] ISO/IEC 29199-2, JPEG XR image coding system
- Part 2: Image coding specification.
[15] I. Daubechies and W. Sweldens, “Factoring wavelet
transforms into lifting steps,” J. Fourier Anal.
Appl., vol.4, no.3, pp.245–267, 1998.
[9] H.J. Meckelburg and D. Lipka, “Fast Hartley transform algorithm,” IET Electron. Lett., vol.21, no.8,
pp.311–313, Apr. 1985.
[16] G. Strang and T. Nguyen, Wavelets and Filter
Banks, Wellesley-Cambridge Press, 1996.
[10] W.H. Chen, C.H. Smith, and S.C. Fralick, “A
fast computational algorithm for the discrete cosine
transform,” IEEE Trans. Commun., vol.25, no.9,
[17] A. Said and W.A. Pearlman, “A new, fast, and efficient image codec based on set partitioning in hierarchical trees,” IEEE Trans. Circuits Syst. Video
pp.1004–1009, Sep. 1977.
Technol., vol.6, no.3, pp.243–250, June 1996.
- 589 -