//+--------------------------------------------------------------------------------------+
//|                                               Fourier Extrapolation of Indicator.mq4 |
//|                                                               Copyright © 2008, gpwr |
//|                                                                   vlad1004@yahoo.com |
//+--------------------------------------------------------------------------------------+
#property copyright "Copyright © 2009, gpwr"
#property indicator_separate_window
#property indicator_buffers 3
#property indicator_color1 Black
#property indicator_color2 Blue
#property indicator_color3 Red
#property indicator_width1 1
#property indicator_width2 2
#property indicator_width3 2

//Global constants
#define pi 3.141592653589793238462643383279502884197169399375105820974944592

//Input parameters
extern int     LastBar  =200;  //Last bar in the past data
extern int     PastBars =500;  //Number of past bars
extern int     FutBars  =200;  //Number of bars to predict 
extern int     HarmNo   =10;   //Number of frequencies; HarmNo=0 computes PastBars harmonics
extern double  FreqTOL  =0.0001;//Tolerance of frequency calculation
//FreqTOL > 0.001 may not converge

//Indicator buffers
double in[];
double pv[];
double fv[];

//Global variables
int np,nf,lb;

int init()
{
   lb=LastBar;
   np=PastBars;
   nf=FutBars;
   if(HarmNo==0 || HarmNo>np) HarmNo=np;
   
   IndicatorBuffers(3);
   SetIndexBuffer(0,in);
   SetIndexBuffer(1,pv);
   SetIndexBuffer(2,fv);
   SetIndexStyle(0,DRAW_LINE,STYLE_SOLID,1);
   SetIndexStyle(1,DRAW_LINE,STYLE_SOLID,2);
   SetIndexStyle(2,DRAW_LINE,STYLE_SOLID,2);
   SetIndexShift(1,-lb);//past data vector 0..np-1; 0 corresponds to bar=lb
   SetIndexShift(2,nf-lb);//future data vector i=0..nf; nf corresponds to bar=lb
   IndicatorShortName("FEoI");  
   return(0);
}

int deinit(){return(0);}

int start()
{
   ArrayInitialize(in,EMPTY_VALUE);
   ArrayInitialize(pv,EMPTY_VALUE);
   ArrayInitialize(fv,EMPTY_VALUE);  
   
//Choose indicator and find the average of its np past values
   double x[]; //stores indicator values
   ArrayResize(x,np);
   double av=0.0;
   for(int i=-lb;i<np;i++)
   {
      in[i+lb]=0.5+iWPR(NULL,0,50,i+lb)/100.0; //change indicator here
      if(i>=0)
      {
         x[i]=in[i+lb];
         av+=x[i];
      }
   }
   av/=np;
   
//Prepare modeled data
   for(i=0;i<np;i++)
   {
      pv[i]=av;
      if(i<=nf) fv[i]=av;
   }

//Fit trigomometric series
   double w,m,c,s;
   for(int harm=1;harm<=HarmNo;harm++)
   {
      Freq(x,w,m,c,s);
      for(i=0;i<np;i++) 
      {
         pv[i]+=m+c*MathCos(w*i)+s*MathSin(w*i);
         if(i<=nf) fv[i]+=m+c*MathCos(w*i)-s*MathSin(w*i);
      }       
   }
   
//Reorder the predicted vector
   for(i=0;i<=(nf-1)/2;i++)
   {
      double tmp=fv[i];
      fv[i]=fv[nf-i];
      fv[nf-i]=tmp;
   } 
   return(0); 
}
//+--------------------------------------------------------------------------------------+
//Quinn and Fernandes algorithm
void Freq(double x[], double& w, double& m, double& c, double& s)
{
   double z[],num,den;
   ArrayResize(z,np);
   double a=0.0;
   double b=2.0;
   z[0]=x[0]-pv[0];
   while(MathAbs(a-b)>FreqTOL)
   {
      a=b;
      z[1]=x[1]-pv[1]+a*z[0];
      num=z[0]*z[1];
      den=z[0]*z[0];
      for(int i=2;i<np;i++)
      {
         z[i]=x[i]-pv[i]+a*z[i-1]-z[i-2];
         num+=z[i-1]*(z[i]+z[i-2]);
         den+=z[i-1]*z[i-1];
      }
      b=num/den;
   }
   w=MathArccos(b/2.0);
   Fit(x,w,m,c,s);
   return;
}
//+-------------------------------------------------------------------------+
void Fit(double x[], double w, double& m, double& c, double& s)
{
   double Sc=0.0;
   double Ss=0.0;
   double Scc=0.0;
   double Sss=0.0;
   double Scs=0.0;
   double Sx=0.0;
   double Sxc=0.0;
   double Sxs=0.0;
   for(int i=0;i<np;i++)
   {
      double cos=MathCos(w*i);
      double sin=MathSin(w*i);
      Sc+=cos;
      Ss+=sin;
      Scc+=cos*cos;
      Sss+=sin*sin;
      Scs+=cos*sin;
      Sx+=(x[i]-pv[i]);
      Sxc+=(x[i]-pv[i])*cos;
      Sxs+=(x[i]-pv[i])*sin;
   }
   Sc/=np;
   Ss/=np;
   Scc/=np;
   Sss/=np;
   Scs/=np;
   Sx/=np;
   Sxc/=np;
   Sxs/=np;
   if(w==0.0)
   {
      m=Sx;
      c=0.0;
      s=0.0;
   }
   else
   {
      //calculating c, s, and m
      double den=MathPow(Scs-Sc*Ss,2)-(Scc-Sc*Sc)*(Sss-Ss*Ss);
      c=((Sxs-Sx*Ss)*(Scs-Sc*Ss)-(Sxc-Sx*Sc)*(Sss-Ss*Ss))/den;
      s=((Sxc-Sx*Sc)*(Scs-Sc*Ss)-(Sxs-Sx*Ss)*(Scc-Sc*Sc))/den;
      m=Sx-c*Sc-s*Ss;
   }
   return;
}