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

//+------------------------------------------------------------------+
//|CGoertzelCycle class for cycle using the Goertzel Algorithm       |
//+------------------------------------------------------------------+
class CGoertzelCycle
  {
private:
   CGoertzel*m_ga;
   double            m_pivot;
   double            m_slope;
   bool              m_detrend;
   bool              m_squaredamp;
   bool              m_flatten;
   uint              m_cycles;
   uint              m_maxper,m_minper;
   double            m_amplitude[];
   double            m_peaks[];
   double            m_cycle[];
   double            m_phase[];
   int               m_peaks_index[];
   double            m_detrended[];
   double            m_flattened[];
   double            m_raw[];


   void              detrend(const double &in_data[], double &out_data[]);
   double            wavepoint(bool use_cycle_strength,uint max_cycles);
   bool              spectrum(const double &in_data[],int shift=-1);
   uint              cyclepeaks(bool use_cycle_strength);

public :
                     CGoertzelCycle(void);
                     CGoertzelCycle(bool detrend,bool squared_amp,bool apply_window, uint min_period,uint max_period);
                    ~CGoertzelCycle(void);


   bool              GetSpectrum(const double &in_data[],double &out_amplitude[]);
   bool              GetSpectrum(bool detrend,bool squared_amp,bool apply_window, uint min_period,uint max_period,const double &in_data[],double &out_amplitude[]);
   uint              GetDominantCycles(bool use_cycle_strength,const double &in_data[],double &out_cycles[]);
   uint              GetDominantCycles(bool use_cycle_strenght,bool detrend,bool squared_amp,bool apply_window, uint min_period,uint max_period,const double &in_data[],double &out_cycles[]);
   void              CalculateDominantCycles(uint prev,uint total,bool use_cycle_strength,const double &in_data[], double &out_signal[]);
   void              CalculateWave(uint prev,uint total,uint max_cycles, bool use_cycle_strength,const double &in_data[], double &out_signal[]);
  };

//+------------------------------------------------------------------+
//|                     Default Constructor                          |
//+------------------------------------------------------------------+
CGoertzelCycle::CGoertzelCycle(void)
  {
   m_pivot=m_slope=0.0;
   m_maxper=m_minper=m_cycles=0;
   m_detrend=m_squaredamp=m_flatten=false;
   m_ga=new CGoertzel();
  }
//+------------------------------------------------------------------+
//|                      Parametric Constructor                      |
//+------------------------------------------------------------------+
CGoertzelCycle::CGoertzelCycle(bool detrend,bool squared_amp,bool apply_window, uint min_period,uint max_period)
  {
   m_pivot=m_slope=0.0;
   m_cycles=0;
   m_maxper=max_period;
   m_minper=min_period;
   m_detrend=detrend;
   m_squaredamp=squared_amp;
   m_flatten=apply_window;
   m_ga=new CGoertzel();

   if(m_maxper)
     {
      // Explicit cast uint -> int for ArrayResize
      ArrayResize(m_amplitude,  (int)m_maxper);
      ArrayResize(m_phase,      (int)m_maxper);
      ArrayResize(m_cycle,      (int)m_maxper);
      ArrayResize(m_peaks,      (int)m_maxper);
      ArrayResize(m_peaks_index,(int)m_maxper);
     }

  }
//+------------------------------------------------------------------+
//|Destructor                                                        |
//+------------------------------------------------------------------+
CGoertzelCycle::~CGoertzelCycle(void)
  {
   if(CheckPointer(m_ga)==POINTER_DYNAMIC)
      delete m_ga;

   ArrayFree(m_amplitude);
   ArrayFree(m_phase);
   ArrayFree(m_cycle);
   ArrayFree(m_peaks_index);
   ArrayFree(m_peaks);
   ArrayFree(m_raw);
   ArrayFree(m_detrended);
   ArrayFree(m_flattened);

  }

//+------------------------------------------------------------------+
//|               helper method for detrending data                  |
//+------------------------------------------------------------------+
void CGoertzelCycle::detrend(const double &in_data[], double &out_data[])
  {
   uint i ;
   double xmean, ymean, x, y, xy, xx ;
   xmean=ymean=x=y=xx=xy=0.0;

   xmean=0.5*double(in_data.Size()-1);

   if(out_data.Size()!=in_data.Size())
      ArrayResize(out_data,in_data.Size()); // in_data.Size() returns int, OK

   for(i=0; i<out_data.Size(); i++)
     {
      x = double(i) - xmean ;
      y = in_data[i] - ymean ;
      xx += x * x ;
      xy += x * y ;
     }

   m_slope=xy/xx;
   m_pivot=xmean;

   for(i=0; i<out_data.Size(); i++)
      out_data[i]=in_data[i]-(double(i)-xmean)*m_slope;
  }
//+-------------------------------------------------------------------------------+
//| helper method that calculates amplitude and phase of spectra of a data sample |
//+-------------------------------------------------------------------------------+
bool CGoertzelCycle::spectrum(const double &in_data[],int shift=-1)
  {
   uint minsize=(3*m_maxper);
   uint fullsize=in_data.Size();
   if(in_data.Size()<minsize)
     {
      Print("Sample size too small in relation to the largest period cycle parameter");
      return false;
     }

   if(!m_ga.SetMinMaxPeriodLength(m_minper,m_maxper))
      return false;

   ArrayInitialize(m_phase,0.0);
   ArrayInitialize(m_amplitude,0.0);
   ArrayInitialize(m_cycle,0.0);
   ArrayInitialize(m_peaks,0.0);
   ArrayInitialize(m_peaks_index,0);

   bool isSeries=ArrayGetAsSeries(in_data);

   uint startshift=(shift<0)?((isSeries)?0:in_data.Size()-1):uint(shift);
   uint startcopy=isSeries?startshift+(minsize-1):startshift-(minsize-1);

   ArrayCopy(m_raw,in_data,0,startcopy,minsize);

   ArraySetAsSeries(m_raw,false);

   if(m_detrend)
      detrend(m_raw,m_detrended);
   else
      ArrayCopy(m_detrended,m_raw);

   double a,b;
   double flat[];

   if(m_flatten)
     {
      if(m_flattened.Size()!=minsize)
         ArrayResize(m_flattened,(int)minsize); // cast uint -> int
      ArrayInitialize(m_flattened,0.0);
      a=m_detrended[0];
      b=(m_detrended[minsize-1]-a)/(minsize-2);
      for(int i=int(minsize-1); i>=0; i--)
         m_flattened[i]=m_detrended[i]-(a+b*(minsize-i));
     }
   else
      ArrayCopy(m_flattened,m_detrended);


   complex dft[];

   if(!m_ga.Dft(m_flattened,dft))
      return false;

   for(uint i=m_minper; i<m_maxper; i++)
     {
      m_amplitude[i]=(m_squaredamp)?(dft[i].real*dft[i].real) + (dft[i].imag*dft[i].imag):MathSqrt((dft[i].real*dft[i].real) + (dft[i].imag*dft[i].imag));
      m_phase[i]=(dft[i].real)?MathArctan(dft[i].imag/dft[i].real):MathArctan(dft[i].imag/1.e-6);

      if(dft[i].real<0)
         m_phase[i]+=M_PI;
      else
        {
         if(dft[i].imag<0)
            m_phase[i]+=2*M_PI;
        }
     }

   m_cycles=0;

   return true;
  }
//+--------------------------------------------------------------------------------+
//| method for identifying spectral peaks and tagging saving them to m_cycles array|
//+--------------------------------------------------------------------------------+
uint CGoertzelCycle::cyclepeaks(bool use_cycle_strength)
  {
   for(uint i=0; i<m_maxper; i++)
     {
      // FIX: explicit cast from uint to int
      m_peaks_index[i] = (int)i;
      if(i>=m_minper+1 && i<m_maxper-1 && m_amplitude[i]>m_amplitude[i+1] && m_amplitude[i]>m_amplitude[i-1])
        {
         m_cycles++;
         m_cycle[i]=double(i);
         m_peaks[i]=use_cycle_strength?m_amplitude[i]:m_amplitude[i]/m_cycle[i];
        }
      else
        {
         m_cycle[i]=0.0;
         m_peaks[i]=-1*DBL_MAX;
        }
     }

   if(m_cycles)
      MathQuickSortDescending(m_peaks,m_peaks_index,0,m_peaks.Size()-1);

   return(m_cycles);
  }
//+------------------------------------------------------------------+
//|method used to calculate filtered wave form using dominant cycles |
//+------------------------------------------------------------------+
double  CGoertzelCycle::wavepoint(bool use_cycle_strength,uint max_cycles)
  {
   double fp,ep;
   fp=ep=0.0;


   if(!m_cycles)
      if(!cyclepeaks(use_cycle_strength))
         return(0);

   uint stop=(max_cycles>m_cycles)?m_cycles:max_cycles;

   uint ss=3*m_amplitude.Size();
   for(uint i = 0; i<stop; i++)
     {
      fp+=m_amplitude[m_peaks_index[i]]*MathCos(2.0*M_PI*MathPow(m_peaks_index[i],-1)*double(ss+1)+m_phase[m_peaks_index[i]]);
      ep+=m_amplitude[m_peaks_index[i]]*MathCos(2.0*M_PI*MathPow(m_peaks_index[i],-1)*double(ss)+m_phase[m_peaks_index[i]]);
     }

   return(fp-ep);

  }
//+------------------------------------------------------------------+
//|          public method to access the amplitude values            |
//+------------------------------------------------------------------+
bool CGoertzelCycle::GetSpectrum(const double &in_data[],double &out_amplitude[])
  {
   if(spectrum(in_data))
     {
      ArrayCopy(out_amplitude,m_amplitude);
      return true;
     }
   return false;
  }
//+------------------------------------------------------------------+
//|         public method to access the amplitude values             |
//+------------------------------------------------------------------+
bool CGoertzelCycle::GetSpectrum(bool detrend,bool squared_amp,bool apply_window,uint min_period,uint max_period,const double &in_data[],double &out_amplitude[])
  {
   m_pivot=m_slope=0.0;
   m_cycles=0;
   m_maxper=max_period;
   m_minper=min_period;
   m_detrend=detrend;
   m_squaredamp=squared_amp;
   m_flatten=apply_window;

   if(m_maxper)
     {
      // Explicit cast uint -> int for ArrayResize
      ArrayResize(m_amplitude,  (int)m_maxper);
      ArrayResize(m_phase,      (int)m_maxper);
      ArrayResize(m_cycle,      (int)m_maxper);
      ArrayResize(m_peaks,      (int)m_maxper);
      ArrayResize(m_peaks_index,(int)m_maxper);
     }

   if(spectrum(in_data))
     {
      ArrayCopy(out_amplitude,m_amplitude);
      return true;
     }
   return false;

  }
//+----------------------------------------------------------------------------+
//|public method to get the dominant cycle periods arranged in descending order|
//+----------------------------------------------------------------------------+
uint CGoertzelCycle::GetDominantCycles(bool use_cycle_strength,const double &in_data[],double &out_cycles[])
  {
   if(!spectrum(in_data))
     {
      ArrayInitialize(out_cycles,0.0);
      return(0);
     }

   cyclepeaks(use_cycle_strength);

   if(out_cycles.Size()!=m_cycles)
      ArrayResize(out_cycles,(int)m_cycles); // cast uint -> int

   for(uint i=0; i<m_cycles; i++)
      out_cycles[i]=m_cycle[m_peaks_index[i]];

   return m_cycles;

  }
//+----------------------------------------------------------------------------+
//|public method to get the dominant cycle periods arranged in descending order|
//+----------------------------------------------------------------------------+
uint CGoertzelCycle::GetDominantCycles(bool use_cycle_strength,bool detrend,bool squared_amp,bool apply_window,uint min_period,uint max_period,const double &in_data[],double &out_cycles[])
  {
   m_pivot=m_slope=0.0;
   m_cycles=0;
   m_maxper=max_period;
   m_minper=min_period;
   m_detrend=detrend;
   m_squaredamp=squared_amp;
   m_flatten=apply_window;

   if(m_maxper)
     {
      // Explicit cast uint -> int for ArrayResize
      ArrayResize(m_amplitude,  (int)m_maxper);
      ArrayResize(m_phase,      (int)m_maxper);
      ArrayResize(m_cycle,      (int)m_maxper);
      ArrayResize(m_peaks,      (int)m_maxper);
      ArrayResize(m_peaks_index,(int)m_maxper);
     }

   return(GetDominantCycles(use_cycle_strength,in_data,out_cycles));

  }
//+-----------------------------------------------------------------------+
//|method used to access calculated dominant cycles , for use in indcators|
//+-----------------------------------------------------------------------+
void CGoertzelCycle::CalculateDominantCycles(uint prev,uint total,bool use_cycle_strength,const double &in_data[], double &out_signal[])
  {

   uint limit =0;
   bool indexOut=ArrayGetAsSeries(out_signal);
   bool indexIn=ArrayGetAsSeries(in_data);

   if(prev<=0)
     {
      uint firstindex=0;
      if(indexOut)
        {
         limit=total-(m_maxper*3);
         firstindex=limit+1;
        }
      else
        {
         limit=m_maxper*3;
         firstindex=limit-1;
        }
      out_signal[firstindex]=0;
     }
   else
     {
      limit=(indexOut)?total-prev:prev;
     }

   uint found_cycles=0;
   if(indexOut)
     {
      for(int ibar=(int)limit; ibar>=0; ibar--)
        {
         spectrum(in_data,(indexIn)?ibar:total-ibar-1);
         out_signal[ibar]=(cyclepeaks(use_cycle_strength))?m_cycle[m_peaks_index[0]]:0.0;
        }
     }
   else
     {
      for(int ibar=(int)limit; ibar<(int)total; ibar++)
        {
         spectrum(in_data,(indexIn)?total-ibar-1:ibar);
         out_signal[ibar]=(cyclepeaks(use_cycle_strength))?m_cycle[m_peaks_index[0]]:0.0;
        }
     }


  }
//+-----------------------------------------------------------------------+
//|method used to access newly formed wave form, for use in indcators     |
//+-----------------------------------------------------------------------+
void CGoertzelCycle::CalculateWave(uint prev,uint total,uint max_cycles,bool use_cycle_strength,const double &in_data[], double &out_signal[])
  {
   uint limit =0;
   bool indexOut=ArrayGetAsSeries(out_signal);
   bool indexIn=ArrayGetAsSeries(in_data);

   if(prev<=0)
     {
      uint firstindex=0;
      if(indexOut)
        {
         limit=total-(m_maxper*3);
         firstindex=limit+1;
        }
      else
        {
         limit=m_maxper*3;
         firstindex=limit-1;
        }
      out_signal[firstindex]=0;
     }
   else
     {
      limit=(indexOut)?total-prev:prev;
     }

   uint found_cycles=0;
   if(indexOut)
     {
      for(int ibar=(int)limit; ibar>=0; ibar--)
        {
         spectrum(in_data,(indexIn)?ibar:total-ibar-1);
         out_signal[ibar]=out_signal[ibar+1]+wavepoint(use_cycle_strength,max_cycles);
        }
     }
   else
     {
      for(int ibar=(int)limit; ibar<(int)total; ibar++)
        {
         spectrum(in_data,(indexIn)?total-ibar-1:ibar);
         out_signal[ibar]=out_signal[ibar-1]+wavepoint(use_cycle_strength,max_cycles);
        }
     }

  }
//+------------------------------------------------------------------+