フーリエ変換(二)

フーリエ変換(二)
――(java)アルゴリズム実装
りさんフーリエへんかん
離散フーリエ変換は数学法をコンピュータ技術と結びつけ,フーリエ変換のような数学ツールの実用化に広い道を開いた.そのため、理論的価値だけでなく、ある意味でより重要な実用価値もあります.
離散フーリエ変換の定義
x(n)がデジタルシーケンスである場合、その離散フーリエ正変換定義は次式で表される
フーリエ逆変換定義は次式で表される
(1)式と(2)式から分かるように,離散フーリエ変換は離散時間信号を直接処理するフーリエ変換である.1つの連続信号をコンピュータデジタル処理する場合は、離散化処理が必要である.これにより,連続信号に対するフーリエ変換の積分過程は自然に求和過程に脱皮する.
高速フーリエ変換(FFT)
高速フーリエ変換は新しい変換ではなく,離散フーリエ変換のアルゴリズムである.この方法は離散フーリエ変換における余分な演算を解析した上で,さらにこれらの繰返し作業を排除する思想指導の下で得られたので,演算において作業量を大幅に節約し,高速演算の目的を達成した.フーリエ変換の使用範囲を拡大した.次に,基本定義から着手し,その原理を議論する.
有限長シーケンス{x(n)}(0<=n
令
したがって、フーリエ変換ペアは、次の式に書くことができる
正変換式(4)を展開すると次のような計算式が得られる.
上の方程式(6)は行列で表すことができる
以上の演算から明らかなように,各周波数成分を得るにはN回乗算とN−1回加算が必要である.変換全体を完了するにはN 2乗算とN(N−1)加算が必要である.シーケンスが長い場合、必然的に多くの時間がかかります.
上記係数行列を見ると、WmnはN周期、すなわち
例えば、N=8の場合、その周期性は図3~6に示す.exp(iθ)=cosθ+isinθ,すなわちexp(-iθ)=cos(-θ)+isin(-θ)=cosθ-isinθ,だから
N=8の場合、
図1,N=8におけるWmnの周期性と対称性
高速フーリエ変換はFFTと略称する.アルゴリズムは分解の特徴によって一般的に2種類あり,1種類は時間分解,1種類は周波数分解である.FFTの基本形式と演算バタフライフローチャートを紹介します.
x(n)を偶数点と奇数点に分ける.
これにより、離散フーリエ変換は、以下のように書くことができる.
(10)
このアルゴリズムのフローチャートは図2に示すように、図(a)入力は順序であり、演算結果は乱順である.図(b)入力は乱順であり、演算結果は順序である.
図2,FFTバタフライ演算フローチャート(時間分解)
コンピュータで高速フーリエ変換を実現
FFTバタフライフローチャートアルゴリズムを用いてコンピュータ上で高速フーリエ変換を実現するには、以下の問題を解決しなければならない.
1)、反復回数rの決定;
バタフライフローチャートから,反復回数rはNに関係することが分かった.値は次の式で決定できます.

int r = (int)(Math.log10(n)/Math.log10(2));	//     r

式中Nは変換シーケンスの長さである.前記基数2に対するバタフライフローチャートは、2の整数次べき乗である.たとえば、シーケンス長が8の場合は3回、シーケンス長が16の場合は4回などです.
2)、対偶ノードの計算;
フローチャートではxl(k)と表記された点をノードと呼ぶ.ここで、下付きlは列数、すなわち数回目の反復であり、例えばxl(k)は、それが最初の反復の結果であることを示す.kは、フローチャートの行数、すなわちシーケンスのシーケンス番号を表す.各ノードの値は、前のノードペアで計算されます.例えば、x 1(0)およびx 1(4)は、いずれもx(0)およびx(4)で計算される.バタフライフローチャートでは,同じソースを持つ一対のノードを対偶ノードと呼ぶ.
例えば、x 1(0)およびx 1(4)は、x(0)およびx(4)に由来する一対の対偶ノードである.対偶ノードの計算は,各反復における対偶ノードの間隔またはピッチを求める.フローチャートから分かるように、1回目の反復のピッチはN/2、2回目の反復のピッチはN/22、3回目の反復のピッチはN/23などである.以上の解析から,次のような対偶ノードの計算方法が得られる.
ノードがxl(k)である場合、その対偶ノードは
xl(k+N/2l) (12)

if(k < n/Math.pow(2, l)) {
	x1 = k;
	x2 = x1 + (int)(n/Math.pow(2, l));
} else {
	x2 = k;
	x1 = x2 - (int)(n/Math.pow(2, l));
}

式中、lは数回目の反復を示す数字であり、kはシーケンスのシーケンス番号数、Nはシーケンス長である.
例:シーケンス長N=8の場合、x 2(1)の対偶ノードを求める.
xl(k+N/2l)=x2((1+8/22)=x2(3)
x2((1)=x1((1)+W80x1(3)
x2((3)=x1((1)+W84x1(3)
3)、重み付け係数WNPの計算;
WNPの計算は主にp値の決定である.p値は以下の方法で求めることができる.
(1)k値をrビットのバイナリ数(kはシーケンスのシーケンス番号数、rは反復回数)と書く.
(2)このバイナリ数をr-lビット右にシフトし、左の空席をゼロに補正する(結果はrビットである).
(3)この右シフト後のバイナリ数をビット反転する.
(4)このビットを逆転させたバイナリ数を10進数にするとp値が得られる.
例:x 2(2)の重み係数W 8 Pを求める.
(1)k=2なのでバイナリ数010と書く.
(2)r-l=3-2=1,010を右にシフトして001を得る.
(3)001をビット順に逆さまにし、すなわちビットを逆さまにして100を得る.
(4)100を10進数に訳して4を得るのでp=4,x 2(2)の重み付け値はWNPとなる.
対偶ノードの計算と併せて、WNPは、あるノードの重み付け係数がWNPであれば、その対偶ノードの重み付け係数は必然的にWNp+N/2であり、しかもWNP=-WNp+N/2であるため、一対の対偶ノードは次式で計算することができる

	/**
	 *      
	 * 1.  k  r      ;2.         r-l ;3. r          ;4.                 ;
	 * @param k         
	 * @param l    
	 * @param r       
	 * @return     
	 */
	private int getWeight(int k, int l, int r) {
		int d = r-l;	//   
		k = k>>d;		
		return reverseRatio(k, r);
	}

4)、並べ替えの問題.
バタフライフローチャートから、シーケンスx(n)が順番に配列されている場合、バタフライ演算後、その変換シーケンスX(m)は非順番に配列され、すなわち乱順であることがわかる.逆に,x(n)が乱順であれば,順序である.したがって、出力の使用を容易にするためには、x(n)とその変換係数X(m)との対応関係を保証するために、再ソートプログラムを加えることが望ましい.具体的なソート方法は次のとおりです.
(1)最後の反復結果xl(k)のシーケンス番号kをバイナリ数とする.
(2)rビットのバイナリビットを逆転する:
(3)逆転後の2進数が表す10進数を求めると,x(k)に対応するX(m)のシーケンス番号数が得られる.
例:
N=8の最終反復結果:
x 3(0)→000→逆転→000→十進法(0)
x 3(1)→001→逆転→100→十進法(4)
x 3(2)→010→逆転→010→十進法(2)
x 3(3)→011→逆転→110→十進法(6)
x 3(4)→100→逆転→001→十進法(1)
x 3(5)→101→逆転→101→十進法(5)
x 3(6)→110→逆転→011→十進法(3)
x 3(7)→111→逆転→111→十進法(7)

	/**
	 *          ，  0101   1010
	 * 1.  k  r      ;2. r          ;3.                 ;
	 * @param k         
	 * @param r       
	 * @return         
	 */
	private int reverseRatio(int k, int r) {
		int n = 0;
		StringBuilder sb = new StringBuilder(Integer.toBinaryString(k));
		StringBuilder sb2 = new StringBuilder("");
		if(sb.length()

一次元フーリエ変換のソースコード実装
複素数の計算が用いられているが,複素数の計算のソースコード実装については,『
シミュレーション複素数とその演算」

/**
	 *          
	 * @param values        
	 * @return           
	 */
	public Complex[] fft(Complex[] values) {
		int n = values.length;
		int r = (int)(Math.log10(n)/Math.log10(2));	//     r
		Complex[][] temp = new Complex[r+1][n];	//         
		Complex w = new Complex(); 	//   
		temp[0] = values;
		int x1, x2;	//          
		int p, t;	//p      Wpn p , t            
		for(int l=1; l<=r; l++) {
			if(l != r) {
				for(int k=0; k

2 Dフーリエ変換のソースコード実装

/**
	 *           
	 * @param matrix        	 
	 * @param w     
	 * @param h     
	* @return           
	 */
	public Complex[][] fft(Complex matrix[][], int w, int h) {
		double r1 = Math.log10(w)/Math.log10(2.0) - (int)(Math.log10(w)/Math.log10(2.0));
		double r2 = Math.log10(h)/Math.log10(2.0) - (int)(Math.log10(w)/Math.log10(2.0));		
		if(r1 != 0.0 || r2 != 0.0) {
			System.err.println("     w h  2 n  ！");
			return null;
		}
		int r = 0;
		r = (int)(Math.log10(w)/Math.log10(2));
		//        
		for(int i=0; i