離散フーリエ変換 – AS3.0

“画像” に対する離散フーリエ変換について。 ”音” は扱わないので注意。
また、一度で書くには分量が多いので何回かに分けて書こうと思います。

・離散フーリエ変換 (DFT)

フーリエ変換についてはFlasherの間でも音関連でよく話題に出ているようなので理論的な説明は省きます。
画像の場合も音の周波数解析と基本は同じ。
ただ、音と違って画像は二次元なので空間周波数領域への変換を行うことになります。
(音の信号 x(t)は一次元、画像の信号 f(x,y) は二次元になる)
実際の計算ではX軸方向にフーリエ変換した後、データを転置してY軸方向にもう一度フーリエ変換します。
また、高速に計算を行うために実装上では高速フーリエ変換(FFT)を用いるのが一般的です。

とりあえずは一次元のFFT(1D-FFT)から始めます。
ActionScriptでの例は以下のようになります。
・1D-FFT

package {
	import __AS3__.vec.Vector;

	public class FFT{
		private var n:int; // データ数
		private var bitrev:Vector.<int>; // ビット反転テーブル
		private var cstb:Vector.<Number>; // 三角関数テーブル

		public function FFT(n:int) {
			if(n != 0 && (n & (n-1)) == 0){
				this.n = n;
				this.cstb = new Vector.<Number>(n+(n>>2), true);
				this.bitrev = new Vector.<int>(n, true);
				makeCstb();
				makeBitrev();
			}  
		}
		// FFT
		public function fft(re:Vector.<Number>, im:Vector.<Number>):Boolean {
			return fftCore(re, im, 1);
		}
		//逆FFT
		public function ifft(re:Vector.<Number>, im:Vector.<Number>):Boolean {
			if(fftCore(re, im, -1)){
				for(var i:int=0; i<n; i++){
					re[i] /= n;
					im[i] /= n;
				}
			}
			return true;
		}
		private function fftCore(re:Vector.<Number>, im:Vector.<Number>, sign:int):Boolean {
			var h:int, d:int, wr:Number, wi:Number, ik:int, xr:Number, xi:Number;
			for(var i:int=0; i<n; i++){ // ビット反転
				var j:int = bitrev[i];
				if(i<j){
					var tmp:Number = re[i];
					re[i] = re[j];
					re[j] = tmp;
					tmp = im[i];
					im[i] = im[j];
					im[j] = tmp;
				}
			}
			for(var k:int=1; k<n; k<<=1){ // バタフライ演算
				h = 0;
				d = n/(k<<1);
				for(var j:int=0; j<k; j++){
					wr = cstb[h+(n>>2)];
					wi = sign*cstb[h];
					for(var i:int=j; i<n; i+=(k<<1)){
						ik = i+k;
						xr = wr*re[ik] + wi*im[ik]
						xi = wr*im[ik] - wi*re[ik];
						re[ik] = re[i] - xr;
						re[i] += xr;
						im[ik] = im[i] - xi;
						im[i] += xi;
					}
					h += d;
				}
			}
			return true;
		}
		private function makeCstb():void {
			var n2:int = n>>1, n4:int = n>>2, n8:int = n>>3;
			var t:Number = Math.sin(Math.PI/n);
			var dc:Number = 2*t*t;
			var ds:Number = Math.sqrt(dc*(2-dc));
			t = 2*dc;
			var c:Number = cstb[n4] = 1;
			var s:Number = cstb[0] = 0;
			for(var i:int=1; i<n8; i++){
				c -= dc;  dc += t*c;
				s += ds;  ds -= t*s;
				cstb[i] = s;  cstb[n4-i] = c;
			}
			if(n8 != 0) cstb[n8] = Math.sqrt(0.5);
			for(var j:int=0; j<n4; j++) cstb[n2-j] = cstb[j];
			for(var k:int=0; k<n2+n4; k++) cstb[k+n2] = -cstb[k];
		}
		private function makeBitrev():void {
			var i:int = 0, j:int = 0, k:int = 0;
			bitrev[0] = 0;
			while(++i<n){
				k = n >> 1;
				while(k<=j){
					j -= k;
					k >>= 1 ;
				}
				j += k;
				bitrev[i] = j;
			}
		}
	}
}

package {

import __AS3__.vec.Vector;

public class FFT{

private var n:int; // データ数

private var bitrev:Vector.<int>; // ビット反転テーブル

private var cstb:Vector.<Number>; // 三角関数テーブル

public function FFT(n:int) {

if(n != 0 && (n & (n-1)) == 0){

this.n = n;

this.cstb = new Vector.<Number>(n+(n>>2), true);

this.bitrev = new Vector.<int>(n, true);

makeCstb();

makeBitrev();

}

// FFT

public function fft(re:Vector.<Number>, im:Vector.<Number>):Boolean {

return fftCore(re, im, 1);

}

//逆FFT

public function ifft(re:Vector.<Number>, im:Vector.<Number>):Boolean {

if(fftCore(re, im, -1)){

for(var i:int=0; i<n; i++){

re[i] /= n;

im[i] /= n;

}

return true;

}

private function fftCore(re:Vector.<Number>, im:Vector.<Number>, sign:int):Boolean {

var h:int, d:int, wr:Number, wi:Number, ik:int, xr:Number, xi:Number;

for(var i:int=0; i<n; i++){ // ビット反転

var j:int = bitrev[i];

if(i<j){

var tmp:Number = re[i];

re[i] = re[j];

re[j] = tmp;

tmp = im[i];

im[i] = im[j];

im[j] = tmp;

}

for(var k:int=1; k<n; k<<=1){ // バタフライ演算

h = 0;

d = n/(k<<1);

for(var j:int=0; j<k; j++){

wr = cstb[h+(n>>2)];

wi = sign*cstb[h];

for(var i:int=j; i<n; i+=(k<<1)){

ik = i+k;

xr = wr*re[ik] + wi*im[ik]

xi = wr*im[ik] - wi*re[ik];

re[ik] = re[i] - xr;

re[i] += xr;

im[ik] = im[i] - xi;

im[i] += xi;

}

h += d;

}

return true;

}

private function makeCstb():void {

var n2:int = n>>1, n4:int = n>>2, n8:int = n>>3;

var t:Number = Math.sin(Math.PI/n);

var dc:Number = 2*t*t;

var ds:Number = Math.sqrt(dc*(2-dc));

t = 2*dc;

var c:Number = cstb[n4] = 1;

var s:Number = cstb[0] = 0;

for(var i:int=1; i<n8; i++){

c -= dc; dc += t*c;

s += ds; ds -= t*s;

cstb[i] = s; cstb[n4-i] = c;

}

if(n8 != 0) cstb[n8] = Math.sqrt(0.5);

for(var j:int=0; j<n4; j++) cstb[n2-j] = cstb[j];

for(var k:int=0; k<n2+n4; k++) cstb[k+n2] = -cstb[k];

}

private function makeBitrev():void {

var i:int = 0, j:int = 0, k:int = 0;

bitrev[0] = 0;

while(++i<n){

k = n >> 1;

while(k<=j){

j -= k;

k >>= 1 ;

}

j += k;

bitrev[i] = j;

}

・テスト用

package {
	import __AS3__.vec.Vector;
	import flash.display.Sprite;
	import flash.utils.getTimer;
	public class Test extends Sprite {	
		private const N:int = 32;
		private var fftTest:FFT;
		public function Test() {
			fftTest = new FFT(N);
			var x1:Vector.<Number> = new Vector.<Number>(N, true);
			var y1:Vector.<Number> = x1.concat();
			var x2:Vector.<Number> = x1.concat();
			var y2:Vector.<Number> = x1.concat();
			var x3:Vector.<Number> = x1.concat();
			var y3:Vector.<Number> = x1.concat();
			for(var i:int=0; i<N; i++){
				x1[i] = x2[i] = 6*Math.cos(6*Math.PI*i/N) + 4*Math.sin(18*Math.PI*i/N);
			}
			var ts:int = getTimer();
			fftTest.fft(x2, y2);  // FFT
			var te:int = getTimer();
			trace("time FFT: "+(te-ts)+" ms");
			for(var i:int=0; i<N; i++){
				x3[i] = x2[i];
				y3[i] = y2[i];
			}
			fftTest.ifft(x3, y3);  // 逆FFT
			trace("データ数: (元データ) (FFT) (逆FFT)");
			for(var i:int=0; i<N; i++){
				trace(i+": ("+x1[i]+", "+y1[i]+") ("+x2[i]+", "+y2[i]+") ("+x3[i]+", " +y3[i]+")");
			}
		}
	}
}

package {

import __AS3__.vec.Vector;

import flash.display.Sprite;

import flash.utils.getTimer;

public class Test extends Sprite {

private const N:int = 32;

private var fftTest:FFT;

public function Test() {

fftTest = new FFT(N);

var x1:Vector.<Number> = new Vector.<Number>(N, true);

var y1:Vector.<Number> = x1.concat();

var x2:Vector.<Number> = x1.concat();

var y2:Vector.<Number> = x1.concat();

var x3:Vector.<Number> = x1.concat();

var y3:Vector.<Number> = x1.concat();

for(var i:int=0; i<N; i++){

x1[i] = x2[i] = 6*Math.cos(6*Math.PI*i/N) + 4*Math.sin(18*Math.PI*i/N);

}

var ts:int = getTimer();

fftTest.fft(x2, y2); // FFT

var te:int = getTimer();

trace("time FFT: "+(te-ts)+" ms");

for(var i:int=0; i<N; i++){

x3[i] = x2[i];

y3[i] = y2[i];

}

fftTest.ifft(x3, y3); // 逆FFT

trace("データ数: (元データ) (FFT) (逆FFT)");

for(var i:int=0; i<N; i++){

trace(i+": ("+x1[i]+", "+y1[i]+") ("+x2[i]+", "+y2[i]+") ("+x3[i]+", " +y3[i]+")");

}

・出力
誤差が出るので逆FFTしたときに全く同じ値には戻りません。
(誤差：10のマイナス14~15乗ほど)

time FFT: 1 ms
データ数: (元データ) (FFT) (逆FFT)
0: (6, 0) (1.3100631690576847e-14, 0) (6, 0)
1: (8.911958795428193, 0) (-3.148595663620103e-15, -1.634119351805764e-15) (8.911958795428195, 2.220446049250315e-16)
2: (0.7653668647301803, 0) (-3.582018523007779e-14, 1.0273603512156518e-14) (0.76536686473018, 5.551115123125783e-17)
3: (-4.496420381306951, 0) (96.00000000000003, -4.440892098500626e-14) (-4.496420381306951, -4.440892098500624e-16)
4: (-1.4142135623730954, 0) (2.074753706171886e-14, 3.110202491371294e-14) (-1.4142135623730956, 0)
5: (-3.6624307503409734, 0) (-1.1638332535152583e-14, 4.013150866834269e-16) (-3.6624307503409748, -2.2204460492503091e-16)
6: (-9.238795325112866, 0) (2.0928516568038727e-14, -1.362553389586755e-14) (-9.238795325112866, 0)
7: (-4.113782686182125, 0) (2.6645352591003757e-15, 2.1316282072803006e-14) (-4.113782686182125, 4.440892098500626e-16)
8: (3.999999999999999, 0) (3.3306690738754696e-15, 3.419486915845482e-14) (3.999999999999999, 0)
9: (2.5530601100531034, 0) (-6.838973831690964e-14, -63.99999999999999) (2.553060110053105, -5.551115123125785e-16)
10: (1.8477590650225717, 0) (-9.893236143066244e-15, -2.623249277195817e-14) (1.8477590650225717, 0)
11: (8.106992614497795, 0) (-9.677949537650423e-15, -1.0256825949718076e-14) (8.106992614497795, -4.440892098500628e-16)
12: (7.071067811865481, 0) (-1.6306644963218234e-14, -2.796183999634539e-14) (7.07106781186548, 0)
13: (-2.155336517113407, 0) (-2.1316282072803006e-14, -2.6645352591003757e-14) (-2.1553365171134082, -7.7715611723761e-16)
14: (-3.826834323650897, 0) (3.3666689002106555e-14, -5.5624060545941617e-14) (-3.826834323650898, 6.661338147750939e-16)
15: (-1.0656765522023512, 0) (5.4662944006227367e-14, -1.584497406700777e-14) (-1.065676552202351, -4.718447854656915e-16)
16: (-5.999999999999996, 0) (2.4646951146678475e-14, 0) (-5.999999999999995, 0)
17: (-8.911958795428195, 0) (5.4662944006227367e-14, 1.584497406700777e-14) (-8.911958795428198, -2.2204460492503114e-16)
18: (-0.7653668647301939, 0) (3.3666689002106555e-14, 5.5624060545941617e-14) (-0.7653668647301936, -5.551115123125783e-17)
19: (4.496420381306947, 0) (-2.1316282072803006e-14, 2.6645352591003757e-14) (4.496420381306947, 4.440892098500628e-16)
20: (1.414213562373091, 0) (-1.6306644963218234e-14, 2.796183999634539e-14) (1.4142135623730911, 0)
21: (3.6624307503409756, 0) (-9.677949537650423e-15, 1.0256825949718076e-14) (3.6624307503409765, 2.220446049250317e-16)
22: (9.23879532511287, 0) (-9.893236143066244e-15, 2.623249277195817e-14) (9.23879532511287, 0)
23: (4.113782686182145, 0) (-6.838973831690964e-14, 63.99999999999999) (4.113782686182146, -4.440892098500626e-16)
24: (-3.999999999999997, 0) (3.3306690738754696e-15, -3.419486915845482e-14) (-3.9999999999999973, 0)
25: (-2.5530601100531056, 0) (2.6645352591003757e-15, -2.1316282072803006e-14) (-2.553060110053107, 5.551115123125781e-16)
26: (-1.847759065022569, 0) (2.0928516568038727e-14, 1.362553389586755e-14) (-1.847759065022569, 0)
27: (-8.106992614497791, 0) (-1.1638332535152583e-14, -4.013150866834269e-16) (-8.106992614497791, 4.440892098500624e-16)
28: (-7.071067811865472, 0) (2.074753706171886e-14, -3.110202491371294e-14) (-7.071067811865471, 0)
29: (2.1553365171133887, 0) (96.00000000000003, 4.440892098500626e-14) (2.1553365171133905, 7.771561172376092e-16)
30: (3.826834323650912, 0) (-3.582018523007779e-14, -1.0273603512156518e-14) (3.826834323650913, -6.661338147750939e-16)
31: (1.065676552202345, 0) (-3.148595663620103e-15, 1.634119351805764e-15) (1.0656765522023448, 4.718447854656915e-16)

time FFT: 1 ms

データ数: (元データ) (FFT) (逆FFT)

0: (6, 0) (1.3100631690576847e-14, 0) (6, 0)

1: (8.911958795428193, 0) (-3.148595663620103e-15, -1.634119351805764e-15) (8.911958795428195, 2.220446049250315e-16)

2: (0.7653668647301803, 0) (-3.582018523007779e-14, 1.0273603512156518e-14) (0.76536686473018, 5.551115123125783e-17)

3: (-4.496420381306951, 0) (96.00000000000003, -4.440892098500626e-14) (-4.496420381306951, -4.440892098500624e-16)

4: (-1.4142135623730954, 0) (2.074753706171886e-14, 3.110202491371294e-14) (-1.4142135623730956, 0)

5: (-3.6624307503409734, 0) (-1.1638332535152583e-14, 4.013150866834269e-16) (-3.6624307503409748, -2.2204460492503091e-16)

6: (-9.238795325112866, 0) (2.0928516568038727e-14, -1.362553389586755e-14) (-9.238795325112866, 0)

7: (-4.113782686182125, 0) (2.6645352591003757e-15, 2.1316282072803006e-14) (-4.113782686182125, 4.440892098500626e-16)

8: (3.999999999999999, 0) (3.3306690738754696e-15, 3.419486915845482e-14) (3.999999999999999, 0)

9: (2.5530601100531034, 0) (-6.838973831690964e-14, -63.99999999999999) (2.553060110053105, -5.551115123125785e-16)

10: (1.8477590650225717, 0) (-9.893236143066244e-15, -2.623249277195817e-14) (1.8477590650225717, 0)

11: (8.106992614497795, 0) (-9.677949537650423e-15, -1.0256825949718076e-14) (8.106992614497795, -4.440892098500628e-16)

12: (7.071067811865481, 0) (-1.6306644963218234e-14, -2.796183999634539e-14) (7.07106781186548, 0)

13: (-2.155336517113407, 0) (-2.1316282072803006e-14, -2.6645352591003757e-14) (-2.1553365171134082, -7.7715611723761e-16)

14: (-3.826834323650897, 0) (3.3666689002106555e-14, -5.5624060545941617e-14) (-3.826834323650898, 6.661338147750939e-16)

15: (-1.0656765522023512, 0) (5.4662944006227367e-14, -1.584497406700777e-14) (-1.065676552202351, -4.718447854656915e-16)

16: (-5.999999999999996, 0) (2.4646951146678475e-14, 0) (-5.999999999999995, 0)

17: (-8.911958795428195, 0) (5.4662944006227367e-14, 1.584497406700777e-14) (-8.911958795428198, -2.2204460492503114e-16)

18: (-0.7653668647301939, 0) (3.3666689002106555e-14, 5.5624060545941617e-14) (-0.7653668647301936, -5.551115123125783e-17)

19: (4.496420381306947, 0) (-2.1316282072803006e-14, 2.6645352591003757e-14) (4.496420381306947, 4.440892098500628e-16)

20: (1.414213562373091, 0) (-1.6306644963218234e-14, 2.796183999634539e-14) (1.4142135623730911, 0)

21: (3.6624307503409756, 0) (-9.677949537650423e-15, 1.0256825949718076e-14) (3.6624307503409765, 2.220446049250317e-16)

22: (9.23879532511287, 0) (-9.893236143066244e-15, 2.623249277195817e-14) (9.23879532511287, 0)

23: (4.113782686182145, 0) (-6.838973831690964e-14, 63.99999999999999) (4.113782686182146, -4.440892098500626e-16)

24: (-3.999999999999997, 0) (3.3306690738754696e-15, -3.419486915845482e-14) (-3.9999999999999973, 0)

25: (-2.5530601100531056, 0) (2.6645352591003757e-15, -2.1316282072803006e-14) (-2.553060110053107, 5.551115123125781e-16)

26: (-1.847759065022569, 0) (2.0928516568038727e-14, 1.362553389586755e-14) (-1.847759065022569, 0)

27: (-8.106992614497791, 0) (-1.1638332535152583e-14, -4.013150866834269e-16) (-8.106992614497791, 4.440892098500624e-16)

28: (-7.071067811865472, 0) (2.074753706171886e-14, -3.110202491371294e-14) (-7.071067811865471, 0)

29: (2.1553365171133887, 0) (96.00000000000003, 4.440892098500626e-14) (2.1553365171133905, 7.771561172376092e-16)

30: (3.826834323650912, 0) (-3.582018523007779e-14, -1.0273603512156518e-14) (3.826834323650913, -6.661338147750939e-16)

31: (1.065676552202345, 0) (-3.148595663620103e-15, 1.634119351805764e-15) (1.0656765522023448, 4.718447854656915e-16)

以前研究で医用画像を扱っていたことがあって、その時に周波数解析の基礎を学びました。
でもその後、画像データの数が確保できずにその研究を断念することに；；
(医用画像はそのまま患者の個人情報になるので取り扱いに注意しなければならない。
そのため医療機関のサポート等がないと、研究で使えるデータを十分確保できないなどの理由)

Flashで医用画像を扱うつもりはもちろんありませんが、基礎だけでも形に残しておこうかと。
画像の周波数解析自体は、医用画像に限らず幅広い分野で用いられている基礎的な処理になります。
今回は一次元のFFTでしたが、次はこれを画像に適用するために二次元のFFTを考えたいと思います。

離散フーリエ変換 – AS3.0

・離散フーリエ変換 (DFT)

あわせて読む:

コメントを残す