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

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

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

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

とりあえずは一次元のFFT(1D-FFT)から始めます。
ActionScriptでの例は以下のようになります。
・1D-FFT
[as]
package {
import __AS3__.vec.Vector;

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

public function FFT(n:int) {
if(n != 0 && (n & (n-1)) == 0){
this.n = n;
this.cstb = new Vector.(n+(n>>2), true);
this.bitrev = new Vector.(n, true);
makeCstb();
makeBitrev();
}
}
// FFT
public function fft(re:Vector., im:Vector.):Boolean {
return fftCore(re, im, 1);
}
//逆FFT
public function ifft(re:Vector., im:Vector.):Boolean {
if(fftCore(re, im, -1)){
for(var i:int=0; i, im:Vector., sign:int):Boolean {
var h:int, d:int, wr:Number, wi:Number, ik:int, xr:Number, xi:Number;
for(var i:int=0; i>2)];
wi = sign*cstb[h];
for(var i:int=j; i>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> 1;
while(k<=j){ j -= k; k >>= 1 ;
}
j += k;
bitrev[i] = j;
}
}
}
}
[/as]
・テスト用
[as]
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. = new Vector.(N, true);
var y1:Vector. = x1.concat();
var x2:Vector. = x1.concat();
var y2:Vector. = x1.concat();
var x3:Vector. = x1.concat();
var y3:Vector. = x1.concat();
for(var i:int=0; 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)

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

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

あわせて読む:

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です