//--------------------------------------------------------
#include  "fft.h"
#include <math.h>
//--------------------------------------------------------
/*
  fft_mode　ビット割り当てによるモード制御
  bit0
  0 内部バッファ（src）をFFTクラスが自己所有とする
  1 内部バッファ（src）は外部所有のポインタとして扱う
*/
FFT::FFT(int Bit,  int fft_mode)
{
  FFT_N = pow(2.0, FFT_BIT=Bit);

  if(!((mode_memo=fft_mode) & 0x01)) {
    src = new Vec2<double>[FFT_N];
  }
  sin_d = new double[FFT_N/2];
  cos_d = new double[FFT_N/2];
  Mask_d = new double[FFT_N];
  init();
}
//--------------------------------------------------------
FFT::~FFT()
{
  if(!(mode_memo & 0x01)) {
    delete  []  src;
  }
  delete  []  sin_d;
  delete  []  cos_d;
  delete  []  Mask_d;
}
//--------------------------------------------------------
void  FFT::init(void)
{
	int	m =	FFT_N/2;
  for(int n=0;  n<m;  n++)  {
    sin_d[n] = sin(2.0*PI*n/FFT_N);    cos_d[n] = cos(2.0*PI*n/FFT_N);
  }

  Mask = MDoNothing;   //何もMASKしないのが規定値
}
//--------------------------------------------------------
//外部の整数データを内部の実数データとして書き込む
void  FFT::SetIntData(int	*p)
{
	for(int n=0;	n<FFT_N;	n++)	{
  	src[n].x = *p++;	src[n].y = 0.0;
  }
}
//--------------------------------------------------------
//内部の複素データの絶対値を外部へ整数データとして読み出す
void  FFT::GetIntData(int	*p)
{
	for(int n=0;	n<FFT_N;	n++)	*p++ = (int)abs(src[n]);
}
//--------------------------------------------------------
//Zero=ゼロｄＢの基準値　Dot10=１０ｄＢあたりのドット数
void  FFT::GetDBCnv(int *p, int Dot10, int Zero)  //dB に変換して得る
{
  double  d, zero;
  zero = 1.0/Zero;
	for(int n=0;	n<FFT_N;	n++, p++)  {
    d=power(src[n]);
  	*p = (int)((d<=0.0001)? -40 : 10.0*log10(d*zero))*Dot10*0.1;
  }
}
//--------------------------------------------------------
/*
  窓の定義によってスペクトルの補正計算をする
  ０からFFT_N/2までが補正され、その上は放置される
  FFTが相関計算のためなどならばこれを使ってはいけない。
*/
void  FFT::SpecComp(void)
{
  double  DC, AC;
  switch(Mask)  {
  case  MDoNothing: DC = FFT_N; AC = FFT_N*0.5; break;
  case  MTriangle:  DC = FFT_N *0.5; AC = FFT_N*16.00/64.0; break;
  case  MSine:      DC = FFT_N *40.73/64.0; AC = FFT_N*20.41/64.0; break;
  case  MHunning:   DC = FFT_N *34.56/64.0; AC = FFT_N*17.28/64.0; break;
  default :   DC = FFT_N *32.00/64.0; AC = FFT_N* 16.00/64.0; break;
  }

  DC = 1.0/DC;  AC = 1.0/AC;
  src[0] *= DC; src[FFT_N/2] *= DC;
	for(int n=1;	n<(FFT_N/2)-1;	n++)	src[n] *= AC;
}
//--------------------------------------------------------
void  FFT::bit_rev(void)
{
  int i,  j,  k;
  double swp;
/* bit reverse */
  j= FFT_N / 2;

  for(i=1 ; i<FFT_N ; i++) {
    k= FFT_N;
    if(j  < i)  {
      swp= src[i].x; src[i].x = src[j].x; src[j].x= swp;
      swp= src[i].y; src[i].y = src[j].y; src[j].y= swp;
    }
    for(k= k/2 ;(j>=k) && (k != 0); j= j-k, k=k/2);
    if(k==0) k=1;
    j= j + k;
  }
}
//--------------------------------------------------------
void  FFT::fft_cnv(int  fdir)
{
  int i,  j,  k,  p,  l,  h,  g,  q;
  double 	a, b, sinp, cosp;

  l= FFT_N;  h= 1;

  for(g= 1 ; g <=  FFT_BIT ; g++)  {
    l= l/2;  k= 0;
    for(q= 1;  q <=  h ; q++)  {
      p= 0;
      for(i= k ; i <=  l+k-1 ; i++)  {
        j= i + l;
        a= src[i].x  - src[j].x;
        b= src[i].y  - src[j].y;
        src[i].x = src[i].x + src[j].x;
        src[i].y = src[i].y + src[j].y;

        if(p== 0)  {
          src[j].x = a;  src[j].y= b;
        } else  {
          cosp= cosf(p);  sinp= (fdir>0)? sinf(p) : -sinf(p);
          src[j].x = a*cosp + b*sinp;
          src[j].y = b*cosp - a*sinp;
        }
        p= p + h;
      }
      k= k + l + l;
    }
    h= h + h;
  }
  bit_rev();
}
//----------------------------------------------------
//現在のマスクをかける
void	FFT::MaskExe(void)
{
	if(Mask != MDoNothing)  {
  	for(int	i=0;	i<FFT_N;	i++)	src[i] *=	Mask_d[i];
  }
}
//----------------------------------------------------
void	FFT::DoNothing(void)
{
	Mask = MDoNothing;
}
//----------------------------------------------------
void	FFT::Triangle(void)
{
	int	i,	j;

	if(Mask != MTriangle)	{
		Mask = MTriangle;
	  for(i=0, j = FFT_N/2; i < FFT_N/2; i++)
    	Mask_d[i+j] = 1.0-(Mask_d[i]	=	(double)i/j);
  }
}
//----------------------------------------------------
void	FFT::Sine(void)
{
	if(Mask != MSine)	{
  	Mask = MSine;
	  for(int i=0; i < FFT_N; i++)	Mask_d[i]	=	sin(PI*(double)i/FFT_N);
  }
}
//----------------------------------------------------
void	FFT::Hunning(void)
{
	if(Mask != MHunning)	{
  	Mask = MHunning;
	  for(int i=0; i < FFT_N; i++)
    	Mask_d[i] = 0.54-0.46*cos(2.0*PI*(double)i/FFT_N);
  }
}
//----------------------------------------------------
void	FFT::Humming(void)
{
	if(Mask != MHumming)	{
  	Mask = MHumming;
	  for(int i=0; i < FFT_N; i++)
    	Mask_d[i]	=	0.5-0.5*cos(2.0*PI*(double)i/FFT_N);
  }
}
//----------------------------------------------------
//FFT→逆FFTしたときの規格化　全部をFFT_Nで除する
void  FFT::RevComp(void)
{
  double  RC = 1.0/FFT_N;
	for(int n=0;	n < FFT_N;	n++)	src[n] *= RC;
}

//--------------------------------------------------------
#ifdef  test4fft
//これ以下は試験用、Mathcadのffttest.mcdと同じものにすること
#include <stdio.h>
#include <stdlib.h>
//--------------------------------------------------------
/*実数領域のみを定義してフーリエ変換する
　信号は区切りの良い連続波としてあたえスペクトルを知るための検査
  窓関数による最終補正も検査する
  DoNothing,Hummingのみ試験した

*/
void  fft_test01(void)
{
  class FFT f32(5);  //32
  V2D v2a[32];

  for(int n=0;  n<32; n++)  {
    v2a[n] = V2D(
      1.0
      +3.0*cos(2.0*PI*4.0*n/32.0)
      +0.5*sin(2.0*PI*7.0*n/32.0)
      +0.75*sin(2.0*PI*16.0*n/32.0)
      +0.25*cos(2.0*PI*16.0*n/32.0)
      , 0.0);
  }

  //DoNothing
  for(int n=0;  n<32; n++)  f32.src[n] = v2a[n];
  WriteCSV("fft_test01_src.txt", f32.src, 32);

  f32.fft_cnv();

  WriteCSV("fft_test01_des.txt", f32.src, 32);

  f32.SpecComp();   //前半分の領域にスペクトルが表現される
  WriteCSV("fft_test01_spec.txt", f32.src, 17);
  //参考　http://www.madlabo.com/mad/edat/mathematic/fft/index.htm
  //真ん中は17個目となる FFT_N/2+1

  //Humming窓
  for(int n=0;  n<32; n++)  f32.src[n] = v2a[n];
  f32.Humming();  f32.MaskExe();
  WriteCSV("fft_test01_specHumming1.txt", f32.src, 32);
  f32.fft_cnv();
  f32.SpecComp();   //前半分の領域にスペクトルが表現される
  WriteCSV("fft_test01_specHumming2.txt", f32.src, 17);
}
//--------------------------------------------------------
void  fft_test02(void)   //FFT、逆FFTの検査
{
  class FFT f32(5);  //32
  V2D v2a[32];

  for(int n=0;  n<32; n++)  {
    v2a[n] = V2D(
      1.0
      +3.0*cos(2.0*PI*4.0*n/32.0)
      +0.5*sin(2.0*PI*7.0*n/32.0)
      +0.75*sin(2.0*PI*16.0*n/32.0)
      +0.25*cos(2.0*PI*16.0*n/32.0)
      , 0.0);
  }
  v2a[10].x = 10; //印をつけておく

  //DoNothing
  for(int n=0;  n<32; n++)  f32.src[n] = v2a[n];
  WriteCSV("fft_test02_src.txt", f32.src, 32);

  f32.fft_cnv();
  f32.fft_cnv(-1);
  f32.RevComp();
  WriteCSV("fft_test02_des.txt", f32.src, 32);


}
//--------------------------------------------------------
//速度検査
  class FFT fsp(16);
  V2D v2a[65536];

void  fft_test03(void)
{
  //Init
  int Fbit = 20;
  int FSiz = 65536;

  for(int n=0;  n<FSiz; n++)  {
    v2a[n] = V2D(
      1.0
      +3.0*cos(2.0*PI*4.0*n/32.0)
      +0.5*sin(2.0*PI*7.0*n/32.0)
      +0.75*sin(2.0*PI*16.0*n/32.0)
      +0.25*cos(2.0*PI*16.0*n/32.0)
      , 0.0);
  }

  for(int m=0;  m<1000; m++)  {
    for(int n=0;  n<FSiz; n++)  fsp.src[n] = v2a[n];
    fsp.fft_cnv();
  }
  
}

//--------------------------------------------------------
//--------------------------------------------------------

#endif  //#ifdef  test4fft
//--------------------------------------------------------