HSFFT.cpp

#include "HSFFT.h"
#include <Windows.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_twfac(int n,int k){
	HSFFT_CMPLEX cmp;
	double buf;
	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_twfac(n,k);

	//ob = ia-ib   (行程1）
	ob->Real  = ia->Real - ib->Real;
	ob->Imaginary = ia->Imaginary - ib->Imaginary;
	
	//ob * eb (行程2)

	//途中計算をobr(実部)、obi(虚部)に代入
	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;
	int cnt;
	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;
}

//FFT・IFFT関数
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);

	//引数に不備が有る場合は0を返す
	if(!num_stage || !input || !output)return 0;

	//FFT・IFFT ステージループ
	for(stage_cnt = 0,index_plus = N; stage_cnt < num_stage ; stage_cnt++){
		
		//バタフライ演算で使用するペアを計算する際に使用するもの
		index_plus >>= 1;	//index_plusを2で割る(1ビット右にシフトで実装)

		//バタフライ演算で使用する指数を計算する際に使用するもの
		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;
}

int HS_rfft_exec(double *input, LPHSFFT_CMPLEX output, unsigned int N){
	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]; 

	//回転因子の符号
	k_x = 1;
	
	//ステージ数の計算
	num_stage  = HS_fft_log2(N);

	//引数に不備が有る場合は0を返す
	if(!num_stage || !input || !output)return 0;

	//FFTステージループ
	for(stage_cnt = 0,index_plus = N; stage_cnt < num_stage ; stage_cnt++){
		
		//バタフライ演算で使用するペアを計算する際に使用するもの
		index_plus >>= 1;	//index_plusを2で割る(1ビット右にシフトで実装)

		//バタフライ演算で使用する指数を計算する際に使用するもの
		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].Real = *(input + idx);
					but_i[1].Real = *(input + idx + index_plus);
					but_i[0].Imaginary = but_i[1].Imaginary = 0;
				}

				//バタフライ演算 実行
				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];
			}
		}
	}

	//並べ替え
	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>
#include <windows.h>
#ifndef __HSFFT_H__
#define __HSFFT_H__

//複素数構造体定義部
typedef struct{
	double Real;
	double Imaginary;
}HSFFT_CMPLEX,*PHSFFT_CMPLEX,*LPHSFFT_CMPLEX;

//プロトタイプ宣言部

//2を底としたlog(対数)関数
unsigned int HS_fft_log2(unsigned int);						

//FFTの回転因子
//twfacは『「Tw」iddle「fac」tor』より
HSFFT_CMPLEX HS_fft_twfac(int,int);								

//FFTのバタフライ計算
int HS_fft_butterfly(LPHSFFT_CMPLEX,LPHSFFT_CMPLEX,int,int);	

//ビットリバース
unsigned __int32 HS_fft_BitReverseU32 (unsigned __int32,int);				

//FFT・IFFT関数
int HS_fft_exec(LPHSFFT_CMPLEX, LPHSFFT_CMPLEX , unsigned int, unsigned int);	
int HS_rfft_exec(double*, LPHSFFT_CMPLEX , unsigned int);



#endif /*__HSFFT_H__*/