FFT関数

HSFFT.cpp

#include "HSFFT.h"
#ifndef __HSFFT_CPP__
#define __HSFFT_CPP__

//2を底としたlog(対数)関数 for FFT・IFFT
unsigned int HS_fft_log2(unsigned int r){
	double log_num  , integer;
	log_num = log((double)r) / log(2.0);
	if(modf(log_num , &integer))return 0;
	else return (unsigned int)integer;
}

//FFTの回転子
HSFFT_CMPLEX	HS_fft_roter(int n,int k){
	HSFFT_CMPLEX cmp;
	double buf;
	cmp.Real = cmp.Imaginary = 0;
	buf = (2.0 * M_PI * k)/n;
	cmp.Real = cos(buf);
	cmp.Imaginary -= sin(buf);
	return cmp;
}

//FFTのバタフライ計算
int	HS_fft_butterfly(LPHSFFT_CMPLEX input,LPHSFFT_CMPLEX output,int n ,int k){
	LPHSFFT_CMPLEX ia , ib , oa , ob;
	HSFFT_CMPLEX eb;
	double obr,obi;
	if(!input || !output)return 0;
	
	//入力用 a,b
	ia = input;
	ib = input + 1;
	//出力用 a,b
	oa = output;
	ob = output + 1;

	/*計算*/
	//oa = ia + ib
	oa->Real = ia->Real + ib->Real;
	oa->Imaginary = ia->Imaginary + ib->Imaginary;

	//回転子
	eb = HS_fft_roter(n,k);

	//ob = ia-ib   (行程1）
	ob->Real  = ia->Real - ib->Real;
	ob->Imaginary = ia->Imaginary - ib->Imaginary;
	
	//ob * eb (行程2)
	obr = eb.Real * ob->Real	;
	obr-= eb.Imaginary * ob->Imaginary	;
	obi = eb.Imaginary * ob->Real	;
	obi+= eb.Real * ob->Imaginary	;
	ob->Real = obr;
	ob->Imaginary = obi;
	return 1;
}

//ビットリバース
unsigned __int32 HS_fft_BitReverseU32 (unsigned __int32 num,int bits){
	unsigned __int32 rev_num=0;
	 int cnt=0;
	if(bits > 32)return 0;
	else if(bits == 0) bits = 32;
	for( rev_num = cnt = 0 ; cnt < bits ; cnt++ )
		if((num & ( 1 << cnt )))
			rev_num += ( 1 <<	 (bits - 1 - cnt ));
	return rev_num;
}


int	HS_fft_exec(LPHSFFT_CMPLEX input, LPHSFFT_CMPLEX output , unsigned int N,unsigned int inverse_flag){
	unsigned int k_x , num_stage , stage_cnt , sub_cnt , item_cnt;
	unsigned int k_plus , index_plus , idx;
	unsigned int index_out;
	HSFFT_CMPLEX but_i[2],but_o[2]; 

	//FFT or IFFT
	k_x = (inverse_flag)? -1 : 1 ;
	
	//ステージ数の計算
	num_stage  = HS_fft_log2(N);
	if(!num_stage || !input || !output)return 0;
	index_plus = N;

	//FFT・IFFT ステージループ
	for(stage_cnt = 0 ; stage_cnt < num_stage ; stage_cnt++){
		index_plus/=2;
		k_plus = 1 << stage_cnt;
		//1ステージ当たりの処理ループ
		for(sub_cnt = 0 ; sub_cnt < k_plus ; sub_cnt++){
			//1アイテム当たりの処理
			for(item_cnt = 0 ; item_cnt  < index_plus ; item_cnt++){

				//インデックス(配列要素番号)
				idx = index_plus*sub_cnt*2 + item_cnt;

				//バタフライ演算・入力
				if(stage_cnt){
					but_i[0] = *(output + idx);
					but_i[1] = *(output + idx + index_plus);
				}else{
					but_i[0] = *(input + idx);
					but_i[1] = *(input + idx + index_plus);
				}

				//バタフライ演算
				HS_fft_butterfly(but_i , but_o , N , k_plus * item_cnt * k_x);

				//バタフライ演算結果
				*(output + idx) = but_o[0];
				*(output + idx + index_plus) = but_o[1];
			}
		}
	}
	
	//IFFTならば正規化をする
	if(inverse_flag){
		for(item_cnt=0 ;item_cnt < N ; item_cnt++){
			(output + item_cnt)->Real  /= N;
			(output + item_cnt)->Imaginary /= N;
		}
	}

	//並べ替え
	for(item_cnt=0 ;item_cnt < N ; item_cnt++){
		index_out = HS_fft_BitReverseU32(item_cnt,num_stage);
		if(index_out > item_cnt){
			but_o[0] = *(output + item_cnt);
			*(output + item_cnt) =  *(output + index_out);
			*(output + index_out) = but_o[0];
		}
	}
	return 1;
}

#endif

HSFFT.h

//インクルード部
#include <stdio.h>
#define _USE_MATH_DEFINES	/* math.h内の定数を使う */
#include <math.h>
#include <memory.h>
#ifndef __HSFFT_H__
#define __HSFFT_H__
//複素数構造体定義部
typedef struct{
	double Real;
	double Imaginary;
}HSFFT_CMPLEX,*PHSFFT_CMPLEX,*LPHSFFT_CMPLEX;

//プロトタイプ宣言部
unsigned int HS_fft_log2(unsigned int);						//2を底としたlog(対数)関数
HSFFT_CMPLEX HS_fft_roter(int,int);								//FFTの回転子
int HS_fft_butterfly(LPHSFFT_CMPLEX,LPHSFFT_CMPLEX,int,int);		//FFTのバタフライ計算

int HS_fft_exec(LPHSFFT_CMPLEX, LPHSFFT_CMPLEX , unsigned int, unsigned int);	//FFT・IFFT関数
unsigned __int32 HS_fft_BitReverseU32 (unsigned __int32,int);				//ビットリバース



#endif /*__HSFFT_H__*/