//+------------------------------------------------------------------+
//|                                                    FFT_class.mqh |
//|                                                           denkir |
//|                                                 denkir@gmail.com |
//+------------------------------------------------------------------+
#property copyright "denkir"
#property link      "denkir@gmail.com"
#include <Complex_class.mqh>
//+------------------------------------------------------------------+
//|                FFT Class definition                              |
//+------------------------------------------------------------------+
class CFFT
  {
public:
   Complex           Input[];  //input array of complex numbers
   Complex           Output[]; //output array of complex numbers
public:
   bool              Forward(const uint N);                                   //direct Fourier transformation
   bool              InverseT(const uint N,const bool Scale=true);            //weighted reverse Fourier transformation
   bool              InverseF(const uint N,const bool Scale=false);           //non-weighted reverse Fourier transformation
   void              setCFFT(Complex &data1[],Complex &data2[],const uint N); //set method (1-st variant)
   void              setCFFT(Complex &data1[],Complex &data2[]);              //set method (2-nd variant)
protected:
   void              Rearrange(const uint N);                                 // regrouping
   void              Perform(const uint N,const bool Inverse);                // implementation of transformation
   void              Scale(const uint N);                                     // weighting
  };
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
void CFFT::Scale(const uint N)
  {
   const double Factor=1./double(N);
//   Scale all data entries
   for(uint Position=0; Position<N;++Position)
      Output[Position].opMultEq(Factor);
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
bool CFFT::InverseT(const uint N,const bool Scale=true)
  {
//   Rearrange
   Rearrange(N);
//   Call FFT implementation
   Perform(N,true);
//   Scale if necessary
   if(Scale)
      Scale(N);
//   Succeeded
   return true;
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
bool CFFT::InverseF(const uint N,const bool Scale=false)
  {
//   Rearrange
   Rearrange(N);
//   Call FFT implementation
   Perform(N,true);
//   Scale if necessary
   if(Scale)
      Scale(N);
//   Succeeded
   return true;
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
void CFFT::setCFFT(Complex &data1[],Complex &data2[])
  {
   const uint N=ArraySize(data1);
   ArrayResize(Input,N);ArrayResize(Output,N);
   for(uint i=0;i<N;i++)
     {
      Input[i].opEqual(data1[i]);
      Output[i].opEqual(data2[i]);
     }
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
void CFFT::setCFFT(Complex &data1[],Complex &data2[],const uint N)
  {
   ArrayResize(Input,N);ArrayResize(Output,N);
   for(uint i=0;i<N;i++)
     {
      Input[i].opEqual(data1[i]);
      Output[i].opEqual(data2[i]);
     }
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
bool CFFT::Forward(const uint N) 
  {
//   Rearrange
   Rearrange(N);
//   Call FFT implementation
   Perform(N,false);
//   Succeeded
   return true;
  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
void CFFT::Rearrange(const uint N)
  {
   uint Target=0; // Data entry position
   for(uint Position=0;Position<N;++Position)//   Process all positions of input signal
     {
      Output[Target].opEqual(Input[Position]);
      uint Mask=N;//   Bit mask
      while(Target &(Mask>>=1))//   While bit is set			
         Target &=~Mask;//   Drop bit		
      Target|=Mask;//   The current bit is 0 - set it
     }

  }
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
void CFFT::Perform(const uint N,const bool Inverse /* = false */)
  {
   const double pi=Inverse ? M_PI : -M_PI;
   for(uint Step=1; Step<N; Step<<=1) //   Iteration through dyads, quadruples, octads and so on...
     {
      const uint Jump=Step<<1;//   Jump to the next entry of the same transform factor
      const double delta=pi/double(Step);//   Angle increment
      const double Sine=sin(delta*0.5);//   Auxiliary sin(delta / 2)
      Complex Multiplier,Factor,factor;
      Multiplier.setComplex(-2.*Sine*Sine,sin(delta));//   Multiplier for trigonometric recurrence
      Factor.setComplex(1.0);//   Start value for transform factor, fi = 0
      for(uint Group=0; Group<Step;++Group)//   Iteration through groups of different transform factor
        {
         for(uint Pair=Group; Pair<N; Pair+=Jump)//   Iteration within group
           {
            const uint Match=Pair+Step;//   Match position
            Complex Product;
            Product.opMult(Factor,Output[Match]);//   Second term of two-point transform
            Output[Match].opMinus(Output[Pair],Product); // Transform for fi + pi
            Output[Pair].opPlusEq(Product); //   Transform for fi
           }
         factor.opEqual(Factor);//   Successive transform factor via trigonometric recurrence
         Factor.opMult(Multiplier,factor);
         Factor.opPlusEq(factor);
        }
     }
  }
//+------------------------------------------------------------------+