といってもこれは簡単に考えつくようなものなのですが、高速化系であり、SIMDを使いやすい形なのでちょっと紹介してみます。
アダマール変換の中でも典型的な方法であるウォルシュ・アダマール変換をやってみます。
このアダマール変換は不可逆の画像圧縮や音楽圧縮などで使う離散コサイン変換と似たような性質を持つ変換であり、
かつ計算が加減算のみで構成されているためハードウェアが簡略化しやすいというメリットがあるため、一部の圧縮形式で利用されています。
実際に自前の圧縮ルーチンにもほとんど使いませんが練習のために入っています。
本来、使用するときにはスクランブル処理が必要になりますが、ここでは扱いません。直接の変換のみです。
アダマール変換については以下の漸化式行列を計算することで得られる計算です。
y1,y2,x1,x2をベクトルとして
H(n)を展開して計算したとき
なお、H(0) = 1として、H(1)はy1,y2,x1,x2をスカラ値として扱うとき
展開して
図を使わずにテキストだけで書こうとすると面倒ですね・・・。texの書き方でも勉強した方がいいでしょうか・・・？
それは置いておいて、こんな感じです。適当な解説ページでも見た方が早いです。
わかるとおり、バタフライ演算の形にすることができます。しかも、この形の場合、要素数は自然に2の累乗になります。
処理の開始
関数の開始は以下のように宣言されているものとします。
//アダマール変換を行う
void Hadamard(float *val,size_t n)
{
size_t i,j,k; __m128 xmm0,xmm1,xmm2; __m128 *addvalue,*subvalue;
__declspec(align(16)) float sign[2][4] = { { 1.0f, -1.0f, 1.0f, -1.0f }, { 1.0f, 1.0f, -1.0f, -1.0f } };
もちろん、valは16byteのアライメントがそろっているものとして扱います。signは符号を合わせるための制御値です。
レベル1の処理
以降は浮動小数として処理していきたいと思います。
レベル1のバタフライはほとんどそのまんまできます。
形としては以下の計算を行います。
y1 = x1 + x2
y2 = x1 – x2
y3 = x3 + x4
y4 = x3 – x4
・・・
これをSSE演算で書いてみると
xmm2 = _mm_load_ps(sign[0]);
for(i = 0,addvalue = (__m128 *)val;i < n;i += 4,addvalue++){
xmm1 = xmm0 = _mm_load_ps(addvalue); //x1 x2 x3 x4
xmm0 = _mm_shuffle_ps(xmm0,xmm1,0xb1); //x1 x2 x3 x4 => x2 x1 x4 x3
xmm1 = _mm_mul_ps(xmm1,xmm2); //x1 x2 x3 x4 => x1 -x2 x3 -x4
xmm0 = _mm_add_ps(xmm0,xmm1);
_mm_store_ps(addvalue,xmm0);
}
それほど難しくないはずです。
レベル2の処理
レベル1とちょっと組み合わせが変わるだけです。
y1 = x1 + x3
y3 = x1 – x3
y2 = x2 + x4
y4 = x2 – x4
・・・
xmm2 = _mm_load_ps(sign[1]);
for(i = 0,addvalue = (__m128i *)val;i < n;i += 4,addvalue++){
xmm1 = xmm0 = _mm_load_ps(addvalue); //x1 x2 x3 x4
xmm0 = _mm_shuffle_ps(xmm0,xmm1,0x4e); //x1 x2 x3 x4 => x3 x4 x1 x2
xmm1 = _mm_mul_ps(xmm1,xmm2); //x1 x2 x3 x4 => x1 x2 -x3 -x4
xmm0 = _mm_add_ps(xmm0,xmm1);
_mm_store_ps(addvalue,xmm0);
}
レベルLの処理
後は任意のレベルの処理を使えます。
m = 2 ^ Lとするとき、p = m*i (0 <= i < n/m)とおいて、
y[p+1] = x[p+1] + x[p+1+m]
y[p+2] = x[p+2] + x[p+2+m]
・・
y[p+1+m] = x[p+1] – x[p+1+m]
y[p+2+m] = x[p+2] – x[p+2+m]
・・
ちょっと形が複雑になりますが、バタフライ的に難しいわけではありません。
実際にこれをコード化してみると、
for(i = 4;i < n;i <<= 1){
for(j = 0;j < n;j += i * 2){
addvalue = (__m128 *)(val + j + 0);
subvalue = (__m128 *)(val + j + i);
for(k = 0;k < i;k += 4,addvalue++,subvalue++){
xmm0 = xmm1 = _mm_load_ps(addvalue); //x[p+1] x[p+2] x[p+3] x[p+4]
xmm2 = _mm_load_ps(subvalue); //x[p+1+m] x[p+2+m] x[p+3+m] x[p+4+m]
xmm0 = _mm_add_ps(xmm0,xmm2);
xmm1 = _mm_sub_ps(xmm1,xmm2);
_mm_store_ps(addvalue,xmm0);
_mm_store_ps(subvalue,xmm1);
}
}
}
SIMDを使ったアダマール変換の高速化
SIMDとバタフライ演算はうまく使うことができるパターンが多いです。ここまできれいなバタフライではないかもしれませんが、
計算の形を見てもらえばわかるとおり、かなりきれいな形です。ほとんど無駄はないです。
これはウォルシュ・アダマール変換がそういう形をしているからなんですが・・・。
ここまで大きなレベルでアダマール変換を行うことはそうそう無いとは思います。参考程度に・・・。