Skip to content

Instantly share code, notes, and snippets.

@hirosof

hirosof/HSFFT.cpp

Created May 16, 2012 11:05
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save hirosof/2709558 to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save hirosof/2709558 to your computer and use it in GitHub Desktop.
Download ZIP
FFT関連関数群
Raw
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=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;
}
//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;
}
//振幅スペクトラムを求める
BOOL HS_fft_To_SwingSpectrum(LPHSFFT_CMPLEX lpfft, double *lpss , unsigned int N){
UINT32 i;
if(!lpss || !lpfft)return 0;
for(i=0;i<N;i++)*(lpss + i) = sqrt(pow((lpfft + i)->Real,2) + pow((lpfft + i)->Imaginary , 2));
return TRUE;
}
//パワースペクトラムを求める
BOOL HS_fft_To_PowerSpectrum(LPHSFFT_CMPLEX lpfft, double *lpss , unsigned int N){
UINT32 i;
if(!lpss || !lpfft)return 0;
for(i=0;i<N;i++)*(lpss + i) = pow((lpfft + i)->Real,2) + pow((lpfft + i)->Imaginary , 2);
return TRUE;
}
UINT32 HS_fft_GetMaxSpectrumIndex(double *lpss , unsigned int N){
UINT32 i,id;
double max_d=0;
//最大値を求める
for(i=0;i<N;i++)if(max_d < *(lpss + i)){
id = i;
max_d = *(lpss + i);
}
return id;
}
//振幅スペクトラムをdB(デシベル)に変換する。
UINT32 HS_SwingSpectrum_To_DeciBel(double *lpss ,double *lpdb , unsigned int N,HSREGUL regulmode){
UINT32 i,id;
double max_d=0;
//最大値を求める
for(i=0;i<N;i++)if(max_d < *(lpss + i)){
id = i;
max_d = *(lpss + i);
}
if(regulmode == HSREGUL_RMS)max_d *=sqrt(2.0);
if(max_d==0){
for(i=0;i<N;i++)*(lpdb + i) = 2;
}else{
for(i=0;i<N;i++){
*(lpdb + i) = *(lpss + i)/max_d;
if(*(lpdb + i))*(lpdb + i) = 20*log10(*(lpdb+i));
else *(lpdb + i) = 2;
}
}
return id;
}
#endif
Raw
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;
//REGULは正規という意味の「Regular」から
//正規化モード
typedef enum{
HSREGUL_MAX=0,
HSREGUL_RMS
}HSREGUL;
//プロトタイプ宣言部
//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);
//振幅スペクトラムを求める
BOOL HS_fft_To_SwingSpectrum(LPHSFFT_CMPLEX,double *,unsigned int);
BOOL HS_fft_To_PowerSpectrum(LPHSFFT_CMPLEX , double *, unsigned int);
//振幅スペクトラムをdB(デシベル)に変換する。
UINT32 HS_SwingSpectrum_To_DeciBel(double *,double *, unsigned int,HSREGUL);
UINT32 HS_fft_GetMaxSpectrumIndex(double * , unsigned int );
#endif /*__HSFFT_H__*/
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment