ウォルシュ・アダマール変換

直交変換として有名なフーリエ変換は、変換核として三角関数を用いる。
一方、ウォルシュ・アダマール変換（Walsh Hadamard Transform: 以下WHT）は、矩形関数を変換核として用いる。
アルゴリズムを最後に述べるが、矩形関数を用いるため、定数倍を除いてWHTは加減算のみとなり、実際はフーリエ変換等より高速に処理を行うことができる。
そのため、一時期は熱心に研究された。
しかし、昨今のハードウェアの進歩により、今日ではそれほど魅力的ではなくなっている。

それでも、幾らかは利用価値はあるし、知識の一つとしても損は無いと思える。
また、WHTはインターネット上で調べても、あまり資料が発見できない。
よって、少しでも必要とするものの助けになれば幸いである。

まず、WHTは次の式で定義される。

この逆変換は、

となる。
ただし、

である。
Nはサンプル点の数である。ｐはNの２のｐ乗の係数であり、ｂi(k)はｋを２値表現にしたときの２の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)