[PR]上記の広告は3ヶ月以上新規記事投稿のないブログに表示されています。新しい記事を書く事で広告が消えます。
直交変換として有名なフーリエ変換は、変換核として三角関数を用いる。
一方、ウォルシュ・アダマール変換(Walsh Hadamard Transform: 以下WHT)は、矩形関数を変換核として用いる。
アルゴリズムを最後に述べるが、矩形関数を用いるため、定数倍を除いてWHTは加減算のみとなり、実際はフーリエ変換等より高速に処理を行うことができる。
そのため、一時期は熱心に研究された。
しかし、昨今のハードウェアの進歩により、今日ではそれほど魅力的ではなくなっている。
それでも、幾らかは利用価値はあるし、知識の一つとしても損は無いと思える。
また、WHTはインターネット上で調べても、あまり資料が発見できない。
よって、少しでも必要とするものの助けになれば幸いである。
まず、WHTは次の式で定義される。
この逆変換は、
となる。
ただし、
である。
Nはサンプル点の数である。pはNの2のp乗の係数であり、bi(k)はkを2値表現にしたときの2のi乗の係数である。
上記の式によって生成した、基底画像の様子を次に示す。
確かにパルス波の重ね合わせであることが分かる。
なお、二次元に拡張した時のWHTとその逆変換は次式のようになる。
高速アルゴリズム
(Fast Walsh Hadamard Transform: FWHT)
N=2のp乗とするとき、まず、
と、二進数表現にする。
次に、i=1, 2,・・・,pとして、
を計算していく。
つまり、次のようなバタフライ演算である。
同様に計算していくと最後は、
とWHTが得られる。
かなり端折ったが、以上がFWHTの概要である。
これをCで書いたものを最後に載せて、このエントリーを〆る。
#include <stdio.h> #include <stdlib.h> #define WHT 0 /* walsh-hadamard transform */ #define IWHT 1 /* inverse walsh-hadamard transform */ void make_bitrev(int n, int bitrev[]) { int i, j, k, n2; n2 = n / 2; i = j = 0; for(;;){ bitrev[i] = j; if(++i >= n) break; k = n2; while(k <= j){ j -= k; k /= 2; } j += k; } } /* 1d fast walsh-hadamard transform function Parameters n :x[]の配列長 x[] :入力信号 sw :switch WHT or IWHT Return value 0 :Normal termination. 1000 :n = 0. 2000 :failed memory allocation. */ int fwht_1d(int n, double x[], int sw) { static int p, last_n = 0, *bitrev = NULL; int i, j, k, u, d, n2, ku; unsigned int mask; double temp; if(n != last_n || n == 0){ last_n = n; if(bitrev != NULL) free(bitrev); if(n == 0) return 1000; bitrev = (int*)malloc(n * sizeof(int)); if(bitrev == NULL){ fprintf(stderr, "error: failed memory allocation.\n"); return 2000; } make_bitrev(n, bitrev); p = 0; for(i = n / 2; i != 0; i /= 2) p++; } d = 2; n2 = n / 2; mask = 0x0001; for(i = 1; i <= p; i++){ for(j = 0; j < n2; j++){ u = j * d + d / 2; for(k = j * d; k < u; k++){ temp = x[k]; ku = k ^ mask; x[k] = temp + x[ku]; x[ku] = temp - x[ku]; } } d *= 2; n2 /= 2; mask = mask << 1; } for(i = 0; i < n; i++){ j = bitrev[i]; if(i < j){ temp = x[i]; x[i] = x[j]; x[j] = temp; } } if(!sw){ for(i = 0; i < n; i++){ x[i] /= n; } } return 0; }
参考文献:
「C言語による[最新]アルゴリズム事典」 奥村 晴彦: 技術評論社(1991)
「高速アルゴリズムと並列信号処理」 谷萩 隆嗣: コロナ社(2000)