//+------------------------------------------------------------------+
//|                                                     Goertzel.mqh |
//|                                  Copyright 2023, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Copyright 2023, MetaQuotes Ltd."
#property link      "https://www.mql5.com"
#include<Math\Stat\Math.mqh>
#include<Graphics\Graphic.mqh>

//+------------------------------------------------------------------+
//| CGoertze class implementing Goertzel algorithm                   |
//+------------------------------------------------------------------+
class CGoertzel
  {
private:
   uint              m_minPeriod;
   uint              m_maxPeriod;
public:
   //constructor
                     CGoertzel(void)
     {
      m_minPeriod=m_maxPeriod=0;
     }
   //destructor
                    ~CGoertzel(void)
     {
     }
   //get methods
   uint              GetMinPeriodLength(void) {   return m_minPeriod;}
   uint              GetMaxPerodLength(void) {   return m_maxPeriod;}
   //set methods
   bool              SetMinMaxPeriodLength(const uint min_wavePeriod, const uint max_wavePeriod);
   // goertzel dft transform methods
   bool              Dft(const double &in_data[],complex &out[]);
   bool              Dft(const double &in_data[], double &out_real[], double &out_imaginary[]);
   bool              Dft(const uint min_wavePeriod,const uint max_wavePeriod,const double &in_data[],complex &out[]);
   bool              Dft(const uint min_wavePeriod,const uint max_wavePeriod,const double &in_data[], double &out_real[], double &out_imaginary[]);
  };

//+------------------------------------------------------------------+
//| Set the Minimum and maximum periods of selected frequency band   |
//+------------------------------------------------------------------+
bool CGoertzel::SetMinMaxPeriodLength(const uint min_wavePeriod,const uint max_wavePeriod)
  {
   if(min_wavePeriod<2 || min_wavePeriod>=max_wavePeriod)
     {
      Print("Critical error min_wavePeriod cannot be less than max_wavePeriod or less than 2");
      return false;
     }

   m_minPeriod=min_wavePeriod;
   m_maxPeriod=max_wavePeriod;

   return true;
  }
//+-----------------------------------------------------------------------------------+
//| Calculate the frequency domain representation of input array as complex array     |
//+-----------------------------------------------------------------------------------+
bool CGoertzel::Dft(const double &in_data[],complex &out[])
  {
   uint minsize=(3*m_maxPeriod);
   uint fullsize=in_data.Size();

   if(fullsize<minsize)
     {
      Print("Sample size too small in relation to the largest period cycle parameter");
      return false;
     }

   if(m_minPeriod>=m_maxPeriod || m_minPeriod<2)
     {
      Print("Critical error: Invalid input parameters :- max_period should be larger than min_period and min_period cannot be less than 2");
      return false;
     }

   if(out.Size()!=m_maxPeriod)
      ArrayResize(out,m_maxPeriod);
   
   double v0,v1,v2,freq,coeff,real,imag;


   for(uint i=0; i<m_maxPeriod; i++)
     {
      if(i<m_minPeriod)
        {
         out[i].imag=out[i].real=0.0;
         continue;
        }

      v0=v1=v2=0.0;
      freq=MathPow(i,-1);
      coeff=2.0*MathCos(2.0*M_PI*freq);
      for(uint k=minsize-1; k>0; k--)
        {
         v0=coeff*v1-v2+in_data[k];
         v2=v1;
         v1=v0;
        }

      real=v1-v2*0.5*coeff;

      imag=v2*MathSin(2*M_PI*freq);

      out[i].real=real;
      out[i].imag=imag;

     }

   return true;
  }
//+------------------------------------------------------------------------------------+
//| Calculate the frequency domain representation of input array as two separate arrays|
//+------------------------------------------------------------------------------------+
bool CGoertzel::Dft(const double &in_data[],double &out_real[],double &out_imaginary[])
  {
   uint minsize=(3*m_maxPeriod);
   uint fullsize=in_data.Size();

   if(fullsize<minsize)
     {
      Print("Sample size too small in relation to the largest period cycle parameter");
      return false;
     }

   if(m_minPeriod>=m_maxPeriod || m_minPeriod<2)
     {
      Print("Critical error: Invalid input parameters :- max_period should be larger than min_period and min_period cannot be less than 2");
      return false;
     }

   if(out_real.Size()!=m_maxPeriod)
      ArrayResize(out_real,m_maxPeriod);

   if(out_imaginary.Size()!=m_maxPeriod)
      ArrayResize(out_imaginary,m_maxPeriod);

   double v0,v1,v2,freq,coeff,real,imag;


   for(uint i=0; i<m_maxPeriod; i++)
     {
      if(i<m_minPeriod)
        {
         out_real[i]=out_imaginary[i]=0.0;
         continue;
        }

      v0=v1=v2=0.0;
      freq=MathPow(i,-1);
      coeff=2.0*MathCos(2.0*M_PI*freq);
      for(uint k=minsize-1; k>0; k--)
        {
         v0=coeff*v1-v2+in_data[k];
         v2=v1;
         v1=v0;
        }

      real=v1-v2*0.5*coeff;

      imag=v2*MathSin(2*M_PI*freq);

      out_real[i]=real;
      out_imaginary[i]=imag;

     }

   return true;
  }

//+-----------------------------------------------------------------------------------+
//| Calculate the frequency domain representation of input array as complex array     |
//+-----------------------------------------------------------------------------------+
bool CGoertzel::Dft(const uint min_wavePeriod,const uint max_wavePeriod,const double &in_data[],complex &out[])
  {
   if(!SetMinMaxPeriodLength(min_wavePeriod,max_wavePeriod))
      return(false);

   return Dft(in_data,out);
  }
//+------------------------------------------------------------------------------------+
//| Calculate the frequency domain representation of input array as two separate arrays|
//+------------------------------------------------------------------------------------+
bool CGoertzel::Dft(const uint min_wavePeriod,const uint max_wavePeriod,const double &in_data[], double &out_real[], double &out_imaginary[])
  {
   if(!SetMinMaxPeriodLength(min_wavePeriod,max_wavePeriod))
      return(false);

   return Dft(in_data,out_real,out_imaginary);

  }