Hatena::ブログ(Diary)

gragraphakt

2009-09-21 [AS3.0] 高速フーリエ変換 ( Fast Fourier Transformation )

高速フーリエ変換 ( Fast Fourier Transsformation )

| 08:41

離散フーリエ変換の対称性に着目して、計算量を大幅に削減することができます。

これを高速フーリエ変換と言います。今回は、高速フーリエ変換のお話...。

のはずだったのですが、きんくまデザイン様が既に AS3.0 用のコードを公開しておられましたので、こちらにほんの少しだけ手を加えます。きんくまデザイン様、ありがとうございます。

きんくまデザイン様のコードは、こちらです。

/**
 * FFT
 * Fast Fourier Transformation 高速フーリエ変換
 * @usage <pre> Fourier.FFT( inputReal, inputImaginary, outputReal, outputImaginary, spectrum );</pre>
 * @param inputReal (Array) 
 * @param inputImaginary (Array)
 * @param outputReal (Array)
 * @param outputImaginary (Array)
 * @param spectrum (Array)
 * @return (Boolean)
 */
public static function FFT( inputReal:Array, inputImaginary:Array, 
			    outputReal:Array, outputImaginary:Array,
			    spectrum:Array ):Boolean {
	var N:uint = inputReal.length;
	var ar:Array, ai:Array;
	var m:int, mh:int, i:int, j:int, k:int, irev:int;
	var wr:Number, wi:Number, xr:Number, xi:Number;
	
	ar = new Array();
	ai = new Array();
	for(i = 0; i < N; i++) {
		ar[i] = inputReal[i];
		ai[i] = inputImaginary[i];
	}
	
	i = 0;
	for (j = 1; j < N - 1; j++) {
		for (k = N >> 1; k > (i ^= k); k >>= 1);
			if (j < i) {
				xr = ar[j];
				xi = ai[j];
				ar[j] = ar[i];
				ai[j] = ai[i];
				ar[i] = xr;
				ai[i] = xi;
			}
	}
	
	var angle:Number = 2 * Math.PI / N;
	for (mh = 1; (m = mh << 1) <= N; mh = m) {
		irev = 0;
		for (i = 0; i < N; i += m) {
			wr = Math.cos(angle * irev);
			wi = Math.sin(angle * irev);
			for (k = N >> 2; k > (irev ^= k); k >>= 1);
			for (j = i; j < mh + i; j++) {
				k = j + mh;
				xr = ar[j] - ar[k];
				xi = ai[j] - ai[k];
				ar[j] += ar[k];
				ai[j] += ai[k];
				ar[k] = wr * xr - wi * xi;
				ai[k] = wr * xi + wi * xr;
			}
		}
	}
	
	for(i = 0; i < N; i++) {
		outputReal[i] = ar[i];
		outputImaginary[i] = ai[i];
		spectrum[i] = Math.sqrt( ar[i]*ar[i] + ai[i]*ai[i]);				
 	}
 	
 	return true;

}

使い方は、以下の通りです。


// 高速フーリエ変換をする、入力(実数と虚数)を配列で用意します。
// 実数の配列の長さは、2の累乗にして下さい。例: 512, 1024, 2048,...
// 虚数が無い場合は、実数と同じ長さで、要素が全て 0 の配列を用意します。

// 高速フーリエ変換後の値を保持する配列を用意します。
var outputReal:Array = new Array();      // 実数部
var outputImaginary:Array = new Array(); // 虚数部
var spectrum:Array = new Array(); // スペクトル

// 実行します。
// FFT メソッドは、Fourier クラスのクラスメソッドとします。
Fourier.FFT( inputReal, inputImaginary, outputReal, outputImaginary, spectrum );

// 結果が、outputReal (実数部)、outputImaginary (虚数部)、 spectrum (スペクトル) に入っています。


以上です。

graphakt.

2009-09-17

離散フーリエ変換( DFT: Discrete Fourier Transformation )

| 15:21

今回は、DFT、離散フーリエ変換とその逆変換。C言語で書かれているサンプルは良く見ますが、AS3.0 で書かれているサンプルはあまり見ないので、作ってみました。数学的な導出は割愛します。

/**
 * DFT
 * Discrete Fourier Transformation 離散フーリエ変換
 * @author graphakt
 * @usage <pre> Fourier.DFT( inputReal, inputImaginary, outputReal, outputImaginary, spectrum ); </pre>
 * @param inputReal (Array) 入力 (実数部)
 * @param inputImaginary (Array) 入力 (虚数部)
 * @param outputReal (Array) 出力(実数部)
 * @param outputImaginary (Array) 出力(虚数部)
 * @param spectrum (Array) スペクトル
 * @return (Boolean) 
 */
public static function DFT( inputReal:Array, inputImaginary:Array,
			    outputReal:Array, outputImaginary:Array,
			    spectrum:Array ):Boolean {		
	var pi:Number = Math.PI;
	if( inputReal.length != inputImaginary.length ) return false;
	var N:int = inputReal.length;
	if( N == 0 ) return false;

 	var real:Number      = 0.0;
	var imaginary:Number = 0.0;
	var angle:Number     = 0.0;

	for(var k:uint = 0; k <= N-1; k++) {
		real = 0.0; imaginary = 0.0;
		for(var n:uint = 0; n <= N-1; n++) {
			angle = - 2 * pi * k * n / N;
			real      += inputReal[n] * Math.cos( angle ) - inputImaginary[n] * Math.sin( angle );
			imaginary += inputReal[n] * Math.sin( angle ) + inputImaginary[n] * Math.cos( angle );
		}
		outputReal[k] = real / N;
		outputImaginary[k] = imaginary / N;
		spectrum[k] = Math.sqrt( real*real + imaginary*imaginary );
 	}

	return true;

}

2009/09/21 改修

使うときは以下のようにします。


// 入力に使う実数の配列と虚数の配列を用意します。
// 虚数が無い場合は 0 で実数と同じ長さを持つ配列を用意します。
// ここでは、実数の配列を inputReal 、虚数の配列を inputImaginary とします。

// 離散フーリエ変換後の結果を保持する配列を準備します。
var outputReal:Array = new Array();       // 実数部
var outputImaginary:Array = new Array();  // 虚数部
var spectrum:Array = new Array();         // スペクトル

// 離散フーリエ変換を実行します。
var Fourier.DFT( inputReal, inputImaginary, outputReal, outputImaginary, spectrum );

// outputReal に実数部、outputImaginary に虚数部、 spectrum にスペクトルが入っています。

逆変換は、上記のプログラムを少し弄るだけで作れます。

/**
 * IDFT
 * Inverse Discrete Fourier Transformation 離散フーリエ逆変換
 * @author graphakt
 * @usage <pre> Fourier.IDFT( inputReal, inputImaginary, outputReal, outputImaginary, spectrum );
 * @param inputReal (Array) 入力信号 (実数部)
 * @param inputImaginary (Array) 入力信号 (虚数部)
 * @param outputReal (Array) 出力(実数部)
 * @param outputImaginary (Array) 出力(虚数部)
 * @param spectrum (Array) スペクトル
 * @return (Boolean) 
 */
public static function IDFT( inputReal:Array, inputImaginary:Array,
			     outputReal:Array, outputImaginary:Array,
		             spectrum:Array ):Boolean {		
	var pi:Number = Math.PI;
	if( inputReal.length != inputImaginary.length ) return false;
	var N:int = inputReal.length;
	if( N == 0 ) return false;

 	var real:Number      = 0.0;
	var imaginary:Number = 0.0;		
 	var angle:Number     = 0.0;

	for(var k:uint = 0; k <= N-1; k++) {
		real = 0.0; imaginary = 0.0;
		for(var n:uint = 0; n <= N-1; n++) {
			angle = 2 * pi * k * n / N;
			real      += inputReal[n] * Math.cos( angle ) - inputImaginary[n] * Math.sin( angle );
			imaginary += inputReal[n] * Math.sin( angle ) + inputImaginary[n] * Math.cos( angle );
		}
		outputReal[k] = real;
		outputImaginary[k] = imaginary;
		spectrum[k] = Math.sqrt( real*real + imaginary*imaginary );
 	}

	return true;
}

2009/09/21 改修

間違いやアドバイス等ございましたら、是非ご指摘下さい。

もし使いたいという奇特な方がいらっしゃいましたら、どうぞ使ってあげて下さい。(コメントを下さると私が喜びます)

graphakt.

[参考文献]

画像数学入門 三角関数フーリエ変換から装置まで, 福田覚, 株式会社東洋書店, 2008年8月25日

2009-09-02 2009/09/02

相関値

| 08:19

「DFT, FFT, DCTの基礎と信号処理応用 やり直しのための信号数学」、三谷政昭、CQ出版株式会社、2008年 9月 1日

p25 - p36 での数式をプログラムに直してみました(プログラムが正しいかどうかはかなり怪しいので、ご了承を^^; )

public function init(): void
{
	// 正規基底ベクトル
	var e0: Array = new Array(  1,  1,  1,  1 );
	var e1: Array = new Array(  1,  1, -1, -1 );
	var e2: Array = new Array(  1, -1, -1,  1 );
	var e3: Array = new Array(  1, -1,  1, -1 );
	// 元の信号
	var original: Array = new Array(12, 2, 0, 2);
	// 正規基底ベクトルで元の信号を線形結合したときの展開係数	
	var F0: Number = innerProduct( original, e0 );
	var F1: Number = innerProduct( original, e1 );
	var F2: Number = innerProduct( original, e2 );
	var F3: Number = innerProduct( original, e3 );
	// g0 : 元の信号を 2 サンプル進めて、信号値を 0.5 倍にしたもの
	var g0: Array = multiple( shift(original.slice(), 2), 0.5 );
	// g1 : 元の信号を 2 サンプル進めて、信号値を 0.5 倍にして、 100 加算したもの
	var g1: Array = add( multiple( shift(original.slice(), 2), 0.5 ), 100 );
	// g2 : 元の信号を 2 サンプル進めて、信号値を 0.5 倍にして、 5 倍にしたもの
	var g2: Array = multiple ( multiple( shift(original.slice(), 2), 0.5 ), 5 );
	// g3 : 元の信号と 2 サンプル進めて、信号値を 0.5 倍にして、 5 倍して、 10 を加算したもの
	var g3: Array = add ( multiple ( multiple( shift(original.slice(), 2), 0.5 ), 5 ), 10 );
	// 元の信号と g0 との相関値
	var R_original_g0: Number = correlation( original, g0 );
	// 元の信号と g1 との相関値
	var R_original_g1: Number = correlation( original, g1 );
	// 元の信号と g2 との相関値 
	var R_original_g2: Number = correlation( original, g2 );
	// 元の信号と g3 との相関値
	var R_original_g3: Number = correlation( original, g3 );
	// __original : 元の信号から、元の信号の平均値を引いた信号値
	var __original: Array = add ( original, -average(original) );
	// __g0 : g0 から、 g0 の平均値を引いた信号値
	var __g0: Array = add( g0, -average(g0) );
	// __g1 : g1 から、 g1 の平均値を引いた信号値
	var __g1: Array = add( g1, -average(g1) );
	// __g2 : g2 から、 g2 の平均値を引いた信号値
	var __g2: Array = add( g2, -average(g2) );
	// __g3 : g3 から、 g3 の平均値を引いた信号値
	var __g3: Array = add( g3, -average(g3) );
	// __original と __g0 との相関値
	var __R_original_g0: Number = correlation( __original, __g0 );
	// __original と __g1 との相関値
	var __R_original_g1: Number = correlation( __original, __g1 );
	// __original と __g2 との相関値
	var __R_original_g2: Number = correlation( __original, __g2 );
	// __original と __g3 との相関値
	var __R_original_g3: Number = correlation( __original, __g3 );
}

private function average( a: Array ) :Number {
	var avg:Number = 0.0;
	for (var i: uint = 0; i < a.length; i++ ) {
		avg += a[i];
	}
	return avg / a.length;
}

private function correlation ( a: Array, b: Array ): Number {
	var value: Number = 0.0;
	value = innerProduct(a, b) / ( norm(a) * norm(b) );
	return value;
}

private function norm( a: Array ) :Number {
	var _norm:Number = 0.0;
	for (var i: uint = 0; i < a.length; i++) {
		_norm += a[i]*a[i];
	}
	return Math.sqrt(_norm / a.length);
}
private function innerProduct( a: Array, b: Array):Number {
	if(a.length != b.length) return NaN;
	var ip:Number = 0.0;
	for (var i: uint = 0; i < a.length; i++) {
		ip += a[i]*b[i];
	}
	return ip/a.length;
}
private function add (_original: Array, a: Number): Array {
	for (var i: uint = 0; i < _original.length; i++ ) {
		_original[i] += a;
	}
	return _original;
}
private function multiple ( _original: Array, a: Number ): Array {
	for ( var i: uint = 0; i < _original.length; i++ ) {
		_original[i] *= a;
	}
	return _original;
}
private function shift( _original: Array, n: uint ): Array {
	for( var i: uint = 0; i < n; i++ ) {
		var shifted: Object = _original.shift();
		_original.push(shifted);
	}
	return _original
}

小数点2桁以下が削られていて、その下が気になり、即席でプログラムを書いて確かめてみた、というだけです。

graphakt.

2009-08-23 論理和(OR)

論理和(OR)

| 07:46

論理和は、2つの条件のうち少なくとも1つが真ならば、真を返します。

つまり、両方とも偽の時だけ、偽が返ります。

trace("true  || false ->->->", true  ||false);
trace("false || true  ->->->", false ||true);
trace("true  || true  ->->->", true  ||true);
trace("false || false ->->->", false ||false);

デバッグ結果はこちらです。

true  || false ->->-> true
false || true  ->->-> true
true  || true  ->->-> true
false || false ->->-> false

■例:hoge || fuga という論理和があった場合

AS3.0 では、まず hoge を Boolean 型に変換します。

hoge の型によって、true / false になるかは異なりますが、基本的に null、NaN、0、""(空のString)以外は、 Boolean の true に変換されると覚えています。(例えば、Number 型は初期化しないと NaN になっています)

var _null = null;        trace(_null || true);
var _NaN = NaN;          trace(_NaN  || true);
var _0 = 0;              trace(_0    || true);
var str:String = "";     trace(str   || true);

デバッグ結果は、こちら。

true
true
true
true

hoge を Boolean に変換した後、 hoge が true かどうかを判定します。もし true の場合は、変換の結果ではなく、変換前を返します。右側の fuga まで処理は行きません。

var hoge:int = 0;
trace(++hoge  || false); // 1

1 は true に変換されるので、判定は true で、後ろの false まで処理は届いていません。

var hoge:String = "hoge";
trace(hoge || false); // hoge

こちらは hoge が出力されます。

fuga も同じです。(fuga が判定されるということは、 hoge が false であったことになります。) まず Boolean 型に変換され、 true / false を判定します。 true の場合は返還前を返します。

var fuga:int = 0;
if(true || ++fuga) trace(fuga); // 0

0 が出力されます。 hoge が true なので、 ++fuga は処理されておらず、初期値の 0 のままです。

var fuga:int = 0;
if(false || ++fuga) trace(fuga); // 1

こちらは 1 が出力されます。 hoge が false なので、 ++fuga が処理され、初期値にプラス1された 1 が Boolean型の true に変換されます。 変換前の 1 が返されます。

■もっと詳しい情報

akihiro kamijo [AS3 の論理積と論理和] http://blogs.adobe.com/akamijo/archives/2008/06/as3_4.html

&&=演算子や||=演算子の解説もあります。



論理和のお話でした。

2009-08-20 浅いコピーと深いコピー (shallow copy and deep copy)Comments

浅いコピーと深いコピー (shallow copy and deep copy)

| 03:26

配列のお話。浅いコピーとか深いコピーについて。

var arr0:Array = new Array();
var arr1:Array = new Array();

arr0.push(1,2,3,4,5);

trace("----------------------------コピー前");
trace("arr0 :", arr0.toString());
trace("arr1 :", arr1.toString());

trace("----------------------------コピー後");				
arr1 = arr0.slice();
trace("arr0 :", arr0.toString());
trace("arr1 :", arr1.toString());

trace("----------------------------コピー後に、配列の[2]番目の要素を操作");								
arr1[2] = 30;
trace("arr0 :", arr0.toString());
trace("arr1 :", arr1.toString()); 

デバッグの結果は、以下の通りです。

----------------------------コピー前
arr0 : 1,2,3,4,5
arr1 : 
----------------------------コピー後
arr0 : 1,2,3,4,5
arr1 : 1,2,3,4,5
----------------------------コピー後に、配列の[2]番目の要素を操作
arr0 : 1,2,3,4,5
arr1 : 1,2,30,4,5  <- コピーした配列の要素を操作しても、コピー元に影響は無い。

オブジェクトである要素が元の配列に格納されていない場合は、コピー後の配列の要素を操作しても、コピー元の配列は影響を受けません。次に、オブジェクトである要素が元の配列に格納されている場合です。

// オブジェクトである要素が元の配列に格納されている場合
var arr0:Array = new Array();
var arr1:Array = new Array();

// Object クラスのインスタンスに hoge というプロパティをくっ付けて、配列に格納。
for(var i:uint=0; i<10; i++) {
	var obj:Object = new Object()
	obj.hoge = i;
	arr0.push(obj);
}

// 配列の中身の確認用。
var _trace:Function = function():void {
	var i:uint;
	var output:String;
	for(i=0,output="";i<arr0.length;i++){ output += arr0[i].hoge.toString() + " "; }
	trace("arr0 : ", output);
	for(i=0,output="";i<arr1.length;i++){ output += arr1[i].hoge.toString() + " "; }
	trace("arr1 : ", output);
}

trace("----------------------------コピー前");
_trace();

arr1 = arr0.slice();
trace("----------------------------コピー後");
_trace();

trace("----------------------------コピー後に、配列の[2]番目の要素を操作");								
arr1[2].hoge = "変更しました";
_trace();

デバッグ結果はこちら。

----------------------------コピー前
arr0 :  0 1 2 3 4 5 6 7 8 9 
arr1 :  
----------------------------コピー後
arr0 :  0 1 2 3 4 5 6 7 8 9 
arr1 :  0 1 2 3 4 5 6 7 8 9 
----------------------------コピー後に、配列の[2]番目の要素を操作
arr0 :  0 1 変更しました 3 4 5 6 7 8 9 
arr1 :  0 1 変更しました 3 4 5 6 7 8 9 <- コピーした配列の要素を操作すると、コピー元に影響がある。

オブジェクトを要素として持っている配列を slice(または concat )でコピーした場合、コピーした配列の要素を操作すると、コピー元にも影響があります。オブジェクトを要素として持っている配列でも、コピー後の配列とコピー元の配列とを別々に扱いたい(お互いの要素を操作しても影響が出ない)ようにするためには、以下のようにします。

private function clone(source:Object):*
{
	var ba:ByteArray = new ByteArray();
	ba.writeObject(source);
	ba.position = 0;
	return(ba.readObject());
}

というメソッドを用意します。このメソッドは、下記のアドレスより、Adobe Flex 3 ヘルプから引用しました。

Adobe Flex 3 ヘルプ [配列クローンの作成] http://livedocs.adobe.com/flex/3_jp/html/help.html?content=10_Lists_of_data_1.html

そして、先ほどのソースコードの中の slice を上記のメソッドに置き換えます。

// arr1 = arr0.slice();	
arr1 = clone(arr0);

デバッグの結果はこちらです。

----------------------------コピー前
arr0 :  0 1 2 3 4 5 6 7 8 9 
arr1 :  
----------------------------コピー後
arr0 :  0 1 2 3 4 5 6 7 8 9 
arr1 :  0 1 2 3 4 5 6 7 8 9 
----------------------------コピー後に、配列の[2]番目の要素を操作
arr0 :  0 1 2 3 4 5 6 7 8 9 
arr1 :  0 1 変更しました 3 4 5 6 7 8 9 <- コピーした配列の要素を操作しても、コピー元に影響がない。

ご覧の通り、コピー後の配列を操作しても、コピー元の配列には影響がありません。

浅いコピーと深いコピーのお話でした。