IDCT_SSE2 - QTheora - Theora Decoder Library Wiki - QTheora - Theora Decoder Library

iDCT の SSE2 実装

Theora では、iDCT の実装方法が規格で定められていて、アルゴリズム的な改良の余地がない。

iDCT の C 言語による実装は、このようなものである。

このコードをよく見てみると、16 bit 符号付き整数の、加算、減算、乗算、が使用さていることが わかる。

SSE2 では、16 bit * 8 での、加算、減算が可能なのでこの機能を使用すれば、8 回 ループする 1D iDCT 演算を一回の実行で行うことができる。(PADDW, PSUBB 命令)

乗算部分では、

 #define MUL(T,X) ((COS[T] * (X)) >> 16)

という処理を行っており、これは 16 bit 定数値に、16 bit 符号付き整数を乗算し、上位 16 bit を 取り出す (65536 で割り算する) という処理である。

SSE2 では、PMULHW という、16 bit 符号付き整数２つを乗算し、上位 16 bit を返すという、 まさにそのままの命令がある。(むしろ、この命令を使用する前提で、アルゴリズムを設計したとみるべきであろう。)

ここで問題になるのが、定数テーブルの、

 static const INT32 COS[8] = {
     65536,
     64277,
     60547,
     54491,
     46341,
     36410,
     25080,
     12785
 };

で、16 bit 符号付き整数の範囲 (-32768 - 32767) に収まらない値が含まれていることである。

32768 以上の数は、２の補数と解釈されるので、符号付きの乗算では、(b >= 32768 の場合)

a * (b - 65536)

と計算されてしまう。

なのだが、この式を変形して考えると、

a * (b - 65536) = a * b - a * 65536

さらに 16 bit 算術右シフト (65536 で除算) を考慮すると、

X = (a * (b - 65536)) / 65536 = (a * b) / 65536 - a
(a * b) / 65536 = X + a

となるので、X を PMULHW 命令の結果と考えると、(a * b) / 65536 を求めるには、 PMULHW の結果に a を加えればよいということが分かる。(b >= 32768 のとき)

以上を踏まえて、1D iDCT 関数をコーディングすると以下のようになった。

 ALIGN(0x10) static const UINT16 COS[8][8] = {
         {     8,     8,     8,     8,     8,     8,     8,     8 }, /* 0 */
         { 64277, 64277, 64277, 64277, 64277, 64277, 64277, 64277 }, /* 1 */
         { 60547, 60547, 60547, 60547, 60547, 60547, 60547, 60547 }, /* 2 */
         { 54491, 54491, 54491, 54491, 54491, 54491, 54491, 54491 }, /* 3 */
         { 46341, 46341, 46341, 46341, 46341, 46341, 46341, 46341 }, /* 4 */
         { 36410, 36410, 36410, 36410, 36410, 36410, 36410, 36410 }, /* 5 */
         { 25080, 25080, 25080, 25080, 25080, 25080, 25080, 25080 }, /* 6 */
         { 12785, 12785, 12785, 12785, 12785, 12785, 12785, 12785 }, /* 7 */
 };
 
 #define MUL1(T,X) _mm_add_epi16(_mm_mulhi_epi16(X, C[T]), X)
 #define MUL0(T,X) _mm_mulhi_epi16(X, C[T])
 
 void IDCT_8_SSE2(
         const INT16* x,
         INT16*       y)
 {
         const __m128i* C = (const __m128i*)COS[0];
         const __m128i* X = (const __m128i*)x;
         __m128i*       Y = (__m128i*)y;
 
         __m128i s0;
         __m128i t0, t1, t2, t3, t4, t5, t6, t7;
 
         /* Stage.1 */
 
         s0 = _mm_add_epi16(X[0], X[4]);
         t0 = MUL1(4, s0);
 
         s0 = _mm_sub_epi16(X[0], X[4]);
         t1 = MUL1(4, s0);
 
         t2 = _mm_sub_epi16(MUL0(6, X[2]), MUL1(2, X[6]));
         t3 = _mm_add_epi16(MUL1(2, X[2]), MUL0(6, X[6]));
 
         t4 = _mm_sub_epi16(MUL0(7, X[1]), MUL1(1, X[7]));
         t5 = _mm_sub_epi16(MUL1(3, X[5]), MUL1(5, X[3]));
 
         t6 = _mm_add_epi16(MUL1(5, X[5]), MUL1(3, X[3]));
         t7 = _mm_add_epi16(MUL1(1, X[1]), MUL0(7, X[7]));
 
         /* Stage.2 */
 
         s0 = _mm_sub_epi16(t4, t5);
         t4 = _mm_add_epi16(t4, t5);
         t5 = MUL1(4, s0);
 
         s0 = _mm_sub_epi16(t7, t6);
         t7 = _mm_add_epi16(t7, t6);
         t6 = MUL1(4, s0);
 
         /* Stage.3 */
 
         s0 = _mm_sub_epi16(t0, t3);
         t0 = _mm_add_epi16(t0, t3);
 
         t3 = _mm_sub_epi16(t1, t2);
         t1 = _mm_add_epi16(t1, t2);
 
         t2 = _mm_sub_epi16(t6, t5);
         t6 = _mm_add_epi16(t6, t5);
 
         /* Stage.4 */
 
         Y[0] = _mm_add_epi16(t0, t7);
         Y[1] = _mm_add_epi16(t1, t6);
         Y[2] = _mm_add_epi16(t3, t2);
         Y[3] = _mm_add_epi16(s0, t4);
         Y[4] = _mm_sub_epi16(s0, t4);
         Y[5] = _mm_sub_epi16(t3, t2);
         Y[6] = _mm_sub_epi16(t1, t6);
         Y[7] = _mm_sub_epi16(t0, t7);
 }

この関数は、垂直方向の 1D iDCT 演算を行う。水平方向の iDCT 演算を行うには、行列の転置処理を行って、垂直方向 演算を行い、係数の位置を戻すためにもう一度、転置処理を行えばよい。

iDCT 演算の前処理である、Dequantize のプロセスで、あらかじめ転置処理を行っておけば、追加の転置処理は１回ですむ。

行列転置処理

8 x 8 の行列転置処理は、以下のようなコードとなる。

アンパック命令である PUNPCKLWD, PUNPCKLDQ, PUNPCKLQDQ 命令、 PUNPCKHWD, PUNPCKHDQ, PUNPCKHQDQ 命令を組み合わせて処理する。

 void Transpose_SSE2(
         const INT16* x,
         INT16*       y)
 {
         const __m128i* X = (const __m128i*)x;
         __m128i*       Y = (__m128i*)y;
 
         __m128i u0, u1, u2, u3, u4, u5, u6, u7;
         __m128i t0, t1, t2, t3, t4, t5, t6, t7;
 
         u0 = _mm_unpacklo_epi16(X[0], X[1]);
         u1 = _mm_unpackhi_epi16(X[0], X[1]);
         u2 = _mm_unpacklo_epi16(X[2], X[3]);
         u3 = _mm_unpackhi_epi16(X[2], X[3]);
         u4 = _mm_unpacklo_epi16(X[4], X[5]);
         u5 = _mm_unpackhi_epi16(X[4], X[5]);
         u6 = _mm_unpacklo_epi16(X[6], X[7]);
         u7 = _mm_unpackhi_epi16(X[6], X[7]);
 
         t0 = _mm_unpacklo_epi32(u0, u2);
         t1 = _mm_unpacklo_epi32(u1, u3);
         t2 = _mm_unpackhi_epi32(u0, u2);
         t3 = _mm_unpackhi_epi32(u1, u3);
         t4 = _mm_unpacklo_epi32(u4, u6);
         t5 = _mm_unpacklo_epi32(u5, u7);
         t6 = _mm_unpackhi_epi32(u4, u6);
         t7 = _mm_unpackhi_epi32(u5, u7);
 
         Y[0] = _mm_unpacklo_epi64(t0, t4);
         Y[1] = _mm_unpackhi_epi64(t0, t4);
         Y[2] = _mm_unpacklo_epi64(t2, t6);
         Y[3] = _mm_unpackhi_epi64(t2, t6);
         Y[4] = _mm_unpacklo_epi64(t1, t5);
         Y[5] = _mm_unpackhi_epi64(t1, t5);
         Y[6] = _mm_unpacklo_epi64(t3, t7);
         Y[7] = _mm_unpackhi_epi64(t3, t7);
 }