FFT(Fast Fourier Transform) 및 IFFT(Inverse Fast Fourier Transform)를 수행하는 클래스

FFT(Fast Fourier Transform) 및 IFFT(Inverse Fast Fourier Transform)를 수행하는 클래스

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// This file is part of MatSeis 1.8
// Copyright (c) 2003 Sandia National Laboratories. All rights reserved.
// 
//     John Merchant, bjmerch@sandia.gov
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

package spis.app.anp.com;

/**
 * FFT(Fast Fourier Transform) 및 IFFT(Inverse Fast Fourier Transform)를 수행하는 클래스 
 */
public class FFT
{
	/** 임시 변수 */
	static int _n = 0;
	/** 임시 변수 */
	static int[] le = new int[0];
	/** 임시 변수 */
	static int[] le1 = new int[0];
	/** 임시 변수 */
	static double[] w_r = new double[0];
	/** 임시 변수 */
	static double[] w_i = new double[0];

	/** 
	 * FFT변환시 수치를 조정하는 메소드.
	 * @param n
	 * 			입력값
	 */
	static private void compTwiddles( int n )
	{
		if ( _n == n )
			return;

		_n = n;

		int ln = (int) ( Math.log( (double)n )/Math.log(2) + 0.5 );

		le = new int[ln+1];
		le1 = new int[ln+1];
		w_r = new double[ln+1];
		w_i = new double[ln+1];

		for (int l = 0; l <= ln; l++)
		{       
			le[l] = (int) (Math.exp( (double)l * Math.log(2) ) + 0.5 );
			le1[l] = le[l] / 2;

			w_r[l] = Math.cos( Math.PI / (double)le1[l] );
			w_i[l] = Math.sin( Math.PI / (double)le1[l] );
		}
	}

	/**
	 * FFT(Fast Fourier Transform)변환을 수행하는 메소드. 원본 보존이 됨.
	 * @param f
	 * 			원본값
	 * @param N
	 * 			FFT변환 길이
	 * @return
	 * 			FFT변환의 결과값. 2차 더블형 배열로 구성되며, 0번째 행은 실수부값을 1번째 행에는 허수부값을 보관함. 
	 */
	public static double[][] fft_1d( double[] f, int N )
	{
		// 2의 배수여부를 체크
		N = NumberSummary.nextPow2(N);
		
		int n = f.length;
		double[][] Z = new double[2][N];
		/* 배열 초기화. */
		/*for (int i = 0; i < Z.length; i++) {
			for (int j = 0; j < Z[i].length; j++) {
				Z[i][j] = 0;
			}
		}*/
		
		for (int i=0; i<n; i++)
			Z[0][i] = f[i];
		
		return fft_1d( Z );
	}

	/**
	 * 실제적으로 FFT변환을 수행하는 메소드. 원본 보존이 되지 않음.
	 * @param array
	 * 			2차 더블형 배열로 구성되며, 0번째 행은 실수부값을 1번째 행에는 허수부값을 보관함.
	 * @return
	 * 			FFT변환의 결과값
	 */
	public static double[][] fft_1d( double[][] array )
	{
		if ( !NumberSummary.isPowerOf2(array[0].length) )
			throw new IllegalArgumentException("Array length must be power of 2");
		
		double  u_r,u_i, t_r,t_i;
		int     ln, nv2, k, l, j, ip, i, n;

		n = array[0].length;

		compTwiddles( n );

		ln = (int)( Math.log( (double)n )/Math.log(2) + 0.5 );
		nv2 = n / 2;
		j = 1;
		for (i = 1; i < n; i++ )
		{
			if (i < j)
			{
				t_r = array[0][i - 1];
				t_i = array[1][i - 1];
				array[0][i - 1] = array[0][j - 1];
				array[1][i - 1] = array[1][j - 1];
				array[0][j - 1] = t_r;
				array[1][j - 1] = t_i;
			}
			k = nv2;
			while (k < j)
			{
				j = j - k;
				k = k / 2;
			}
			j = j + k;
		}

		for (l = 1; l <= ln; l++) /* loops thru stages */
		{
			u_r = 1.0;
			u_i = 0.0;
			for (j = 1; j <= le1[l]; j++) /* loops thru 1/2 twiddle values per stage */
			{
				for (i = j; i <= n; i += le[l]) /* loops thru points per 1/2 twiddle */
				{
					ip = i + le1[l];
					t_r = array[0][ip - 1] * u_r - u_i * array[1][ip - 1];
					t_i = array[1][ip - 1] * u_r + u_i * array[0][ip - 1];

					array[0][ip - 1] = array[0][i - 1] - t_r;
					array[1][ip - 1] = array[1][i - 1] - t_i; 

					array[0][i - 1] =  array[0][i - 1] + t_r;
					array[1][i - 1] =  array[1][i - 1] + t_i;  
				} 
				t_r = u_r * w_r[l] + w_i[l] * u_i;
				u_i = w_r[l] * u_i - w_i[l] * u_r;
				u_r = t_r;
			} 
		}  
		return array;
	} /* end of FFT_1d */

	/**
	 * FFT역변환을 수행함. 원본 보존이 되지 않음.
	 * @param array
	 * 			2차 더블형 배열 구조 이며, 0번째 행은 실수부값을 1번째 행에는 허수부값을 보관함.
	 * @return
	 * 			FFT역변환의 결과.
	 */
	public static double[][] ifft_1d( double[][] array )
	{
		double  u_r,u_i, t_r,t_i;
		int     ln, nv2, k, l, j, ip, i, n;

		n = array[0].length;
		compTwiddles(n);

		ln = (int)( Math.log( (double)n )/Math.log(2) + 0.5 );
		nv2 = n / 2;
		j = 1;
		for (i = 1; i < n; i++ )
		{
			if (i < j)
			{
				t_r = array[0][i - 1];
				t_i = array[1][i - 1];
				array[0][i - 1] = array[0][j - 1];
				array[1][i - 1] = array[1][j - 1];
				array[0][j - 1] = t_r;
				array[1][j - 1] = t_i;
			}
			k = nv2;
			while (k < j)
			{
				j = j - k;
				k = k / 2;
			}
			j = j + k;
		}

		for (l = 1; l <= ln; l++) /* loops thru stages */
		{       
			u_r = 1.0;
			u_i = 0.0;
			for (j = 1; j <= le1[l]; j++) /* loops thru 1/2 twiddle values per stage */
			{
				for (i = j; i <= n; i += le[l]) /* loops thru points per 1/2 twiddle */
				{
					ip = i + le1[l];
					t_r = array[0][ip - 1] * u_r - u_i * array[1][ip - 1];
					t_i = array[1][ip - 1] * u_r + u_i * array[0][ip - 1];

					array[0][ip - 1] = array[0][i - 1] - t_r;
					array[1][ip - 1] = array[1][i - 1] - t_i; 

					array[0][i - 1] =  array[0][i - 1] + t_r;
					array[1][i - 1] =  array[1][i - 1] + t_i;  
				} 
				t_r = u_r * w_r[l] - w_i[l] * u_i;
				u_i = w_r[l] * u_i + w_i[l] * u_r;
				u_r = t_r;
			} 
		}  
		return array;
	} /* end of ifft_1d */
	
	/**
	 * FFT역변환을 수행함. 다만, 역변환시 FFT변환전의 원본값과 일치할수 있도록 값을 보정함. 원본 보존이 되지 않음.
	 * @param array
	 * 			2차 더블형 배열 구조 이며, 0번째 행은 실수부값을 1번째 행에는 허수부값을 보관함.
	 * @return
	 * 			FFT역변환의 결과.
	 */
	public static double[][] ifft_1d_n( double[][] array )
	{
		ifft_1d(array);
		
		int length = array[0].length;
		
		// Normalize
		for (int i = 0; i < length; i++) {
			array[0][i] = array[0][i] / (double)length;
			array[1][i] = 0;
		}
		
		return array;
	}
}