#property copyright "www,forex-tsd.com"
#property link      "www,forex-tsd.com"

#property indicator_chart_window
#property indicator_buffers   3
#property indicator_color1    DarkGray
#property indicator_width1    1
#property indicator_style1    3
#property indicator_color2    LimeGreen
#property indicator_width2    2
#property indicator_color3    Orange
#property indicator_width3    2

//---- external variables
extern int     paFilter       =   1;
extern double  paCycles       = 1.0;
extern int     Price          =   0;   //Apply to Price(0-Close;1-Open;2-High;3-Low;4-Median price;5-Typical price;6-Weighted Close)
extern int     Degree         =   1;   //Degree of a Polynomial(no more 12)
extern int     PreSmooth      =   1;   //Period of Preliminary Smoothing
extern int     PreSmoothMode  =   0;   //Mode of MA:0-SMA,1-EMA,2-Wilders(SMMA),3-LWMA
extern double  K_Sigma        =   1;   //Multiplier of Sigma(Standard Deviation)
extern int     FitMode        =   1;   //Fitting Mode: 0-Fitting,1-Moving
extern int     CenterMode     =   1;   //Center Line Plot: 0-without,1-with 
extern int     BarsBack       =   0;   //Shift in history 
extern int     CountBars      =   0;   //Maximum Number of Bars(0-all Bars) 
//---- indicator buffers
double Polynom[];
double UpBand[];
double DnBand[];


//---- global variables
double PriceArray[];
int cBars;
datetime pUpTime = 0, pDnTime = 0, pTime;
int outbars=0, inbars=0, onbars=0;
bool fTime = true;
int Length=20;
double lengths[];
//+---------------------------------------------------------------------------+
//| Custom indicator initialization function                                  |
//+---------------------------------------------------------------------------+
int init()
{
   IndicatorBuffers(4);
   if(BarsBack < 0) BarsBack = 0;
   //IndicatorDigits(MarketInfo(Symbol(), MODE_DIGITS )+2);
   
   SetIndexStyle(0, DRAW_LINE);
   SetIndexStyle(1, DRAW_LINE);
   SetIndexStyle(2, DRAW_LINE);
      
   SetIndexBuffer(0, Polynom);
   SetIndexBuffer(1, UpBand);
   SetIndexBuffer(2, DnBand);
   SetIndexBuffer(3, lengths); SetIndexEmptyValue(3,0);
      
   if(CountBars == 0) cBars = Bars - (Length+PreSmooth-2);
   else cBars = CountBars+1;
   
   SetIndexDrawBegin( 0,Bars-cBars-2);//+2*(Length+PreSmooth));
   SetIndexDrawBegin( 1,Bars-cBars-2);//+2*(Length+PreSmooth));
   SetIndexDrawBegin( 2,Bars-cBars-2);//+2*(Length+PreSmooth));
     
   string short_name="PolyFitBands("+Price+","+Length+","+Degree+")";
   IndicatorShortName(short_name);
   SetIndexLabel(0,"PolyFitCenter");
   SetIndexLabel(1,"UpPolyBand");
   SetIndexLabel(2,"DnPolyBand");
   
   ArrayResize(PriceArray, Length);
      
   return(0);
}

double PolyFit(int mode,double& price[],int deg,int len,int bar)
{
double   result=0,Sum=0;
double   AX[12,12],BX[12],ZX[12],Pow[12];
int      j,k,Row[12];  
static double CX[12];
      
   if (len<=1) Sum = price[bar];
   else
   {	
      if(mode == 1 || (mode == 0 && bar == -len+1))
      {
      for (j=1;j<=deg+1;j++) BX[j]=0;
		
         for (k=1;k<=len;k++) 
  	      {
  	      double YK=price[len-k];
  	      double XK=k;
  	      double Prod=1;
  	         for (j=1;j<=deg+1;j++)
            {
    	      BX[j]+= YK*Prod;
    	      Prod*=XK;
  		      }
         }
      
      for (j=0;j<=2*deg;j++) Pow[j]=0;
	
	   Pow[0]=len;
	
         for (k=1;k<=len;k++)
	      {
  	      XK=k;
  	      Prod=k;
  	         for (j=1;j<=2*deg;j++)
  	         {
    	      Pow[j]+= Prod;    	   
    	      Prod*= XK;
  		      }
         }	
	
	      for (j=1;j<=deg+1;j++)
	      {	
            for (int l=1;l<=deg+1;l++) 
  	         AX[j,l]=Pow[j+l-2];
         }	

	  for (j=1;j<=deg+1;j++) Row[j]=j;
  	
         for (int i=1;i<=deg;i++)
  	      {
  	         for (k=i+1;k<=deg+1;k++)
            {
               if(MathAbs(AX[Row[k],i]) > MathAbs(AX[Row[i],i]))
      	      {	
      	      int Temp=Row[i];
      	      Row[i]=Row[k];
      	      Row[k]=Temp;
               }
            }

            for (k=i+1;k<=deg+1;k++)
    	      {
    	      if(AX[Row[i],i]!=0) AX[Row[k],i]=AX[Row[k],i]/AX[Row[i],i];
    		      for (l=i+1;l<=deg+1;l++)
      	      AX[Row[k],l]=AX[Row[k],l]-AX[Row[k],i]*AX[Row[i],l];
    		   }
         }

      ZX[1]=BX[Row[1]];
      
         for (k=2;k<=deg+1;k++)
  	      {
  	      Sum=0;
         for (l=1;l<=k-1;l++) Sum+=AX[Row[k],l]*ZX[l];
  	   
  	      ZX[k]=BX[Row[k]]-Sum;
	      }
     
      if(AX[Row[deg+1],deg+1] != 0) CX[deg+1]=ZX[deg+1]/AX[Row[deg+1],deg+1];
      	
         for (k=deg;k>=1;k--)
  	      {
  	      Sum=0;
	    		
  	      for (l=k+1;l<=deg+1;l++) Sum += AX[Row[k],l]*CX[l]; 
  	      CX[k]=(ZX[k]-Sum)/AX[Row[k],k];
	      }
      }
      	  
	Sum=CX[deg+1];
	for (k=deg;k>=1;k--) Sum=CX[k]+Sum*(len + bar);
   }
   return (Sum);
}	
//+------------------------------------------------------------------+
//| PolyFitBands_v2.3                                                     |
//+------------------------------------------------------------------+
int start()
{
   int i, j, shift, counted_bars=IndicatorCounted(), limit;
        
   
   if ( counted_bars < 0 )  return(0);
   if ( counted_bars ==0 )  limit=cBars-1; 
         
   if (FitMode==1 && counted_bars > 0) {limit = Bars - counted_bars + BarsBack; int len = Length+PreSmooth;}
   else
   if (FitMode==0) {limit = Length + BarsBack; len = cBars - Length+1;} 
   
   if ( counted_bars < 1 || FitMode==0)
   { 
      for(i=1;i<len;i++)
      { 
      Polynom[cBars-i]=EMPTY_VALUE; 
      UpBand[cBars-i]=EMPTY_VALUE;
      DnBand[cBars-i]=EMPTY_VALUE;
      }
   }   

   if(counted_bars<0) return(-1);
   if(counted_bars>0) counted_bars--;
       int fill = Bars - counted_bars;
      
   for(shift=fill;shift>=0;shift--)
      lengths[shift] = iHilbertPhase(iMA(NULL,0,1,0,MODE_SMA,Price,shift),paFilter,paCycles,shift,0);
   for(shift=BarsBack;shift<limit;shift++ )
   {
      Length = lengths[shift];
         if (ArraySize(PriceArray)!=Length) ArrayResize(PriceArray, Length);
  
      //
      //
      //
      //
      //
      
      if(shift < cBars)
      {
         if((shift==BarsBack && FitMode==0) || (FitMode==1))
         {   
            for(j=0;j<Length;j++) 
            PriceArray[j] = iMA(NULL,0,PreSmooth,0,PreSmoothMode,Price,shift*FitMode + j + (1-FitMode)*BarsBack);
      
            double Sum=0;
            for(j=Length-1;j>=0;j--)
            {
            double poly = PolyFit(0,PriceArray,Degree,Length,-j); 
            double del = PriceArray[j] - poly; 
            Sum += del*del;
            if (FitMode==0) Polynom[j+BarsBack] = poly;  
            }
      
         if(FitMode==1) Polynom[shift] = poly; 
      
            if (Length-1 > 0 && Sum > 0)
            {
               if(Length < 32) double StdDev = MathSqrt(Sum/(Length-1));   
               else
               StdDev = MathSqrt(Sum/Length);       
            }
         }
      
      UpBand[shift] = Polynom[shift] + K_Sigma * StdDev;    
      DnBand[shift] = Polynom[shift] - K_Sigma * StdDev;
      if(CenterMode==0) Polynom[shift] = EMPTY_VALUE;
      }
   }             
   return(0);
}
//+---------------------------------------------------------------------------+




//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
//
//
//
//
//

double workHil[][9];
#define _price      0
#define _smooth     1
#define _detrender  2
#define _period     3
#define _instPeriod 4
#define _phase      5
#define _deltaPhase 6
#define _Q1         7
#define _I1         8

#define Pi 3.14159265358979323846264338327950288

//
//
//
//
//

double iHilbertPhase(double price, double filter, double cyclesToReach, int i, int s=0)
{
   if (ArrayRange(workHil,0)!=Bars) ArrayResize(workHil,Bars);
   int r = Bars-i-1; s = s*9;
      workHil[r][s+_price] = price;
      if (r<50) return(1);
      
   //
   //
   //
   //
   //
      
      workHil[r][s+_smooth]     = (4.0*workHil[r][s+_price]+3.0*workHil[r-1][s+_price]+2.0*workHil[r-2][s+_price]+workHil[r-3][s+_price])/10.0;
      workHil[r][s+_detrender]  = calcComp(r,_smooth,s);
      workHil[r][s+_Q1]         = 0.15*calcComp(r,_detrender,s)  +0.85*workHil[r-1][s+_Q1];
      workHil[r][s+_I1]         = 0.15*workHil[r-3][s+_detrender]+0.85*workHil[r-1][s+_I1];
      workHil[r][s+_phase]      = workHil[r-1][s+_phase];
      workHil[r][s+_instPeriod] = workHil[r-1][s+_instPeriod];

      //
      //
      //
      //
      //
           
         if (MathAbs(workHil[r][s+_I1])>0)
                     workHil[r][s+_phase] = 180.0/Pi*MathArctan(MathAbs(workHil[r][s+_Q1]/workHil[r][s+_I1]));
           
         if (workHil[r][s+_I1]<0 && workHil[r][s+_Q1]>0) workHil[r][s+_phase] = 180-workHil[r][s+_phase];
         if (workHil[r][s+_I1]<0 && workHil[r][s+_Q1]<0) workHil[r][s+_phase] = 180+workHil[r][s+_phase];
         if (workHil[r][s+_I1]>0 && workHil[r][s+_Q1]<0) workHil[r][s+_phase] = 360-workHil[r][s+_phase];

      //
      //
      //
      //
      //
                        
      workHil[r][s+_deltaPhase] = workHil[r-1][s+_phase]-workHil[r][s+_phase];

         if (workHil[r-1][s+_phase]<90 && workHil[r][s+_phase]>270)
             workHil[r][s+_deltaPhase] = 360+workHil[r-1][s+_phase]-workHil[r][s+_phase];
             workHil[r][s+_deltaPhase] = MathMax(MathMin(workHil[r][s+_deltaPhase],60),7);
      
            //
            //
            //
            //
            //
                  
            double alpha    = 2.0/(1.0+MathMax(filter,1));
            double phaseSum = 0; for (int k=0; phaseSum<cyclesToReach*360 && (r-k)>0; k++) phaseSum += workHil[r-k][s+_deltaPhase];
         
               if (k>0) workHil[r][s+_instPeriod]= k;
                  workHil[r][s+_period] = workHil[r-1][s+_period]+alpha*(workHil[r][s+_instPeriod]-workHil[r-1][s+_period]);
   return (workHil[r][s+_period]);
}

//
//
//
//
//

double calcComp(int r, int from, int s)
{
   return((0.0962*workHil[r  ][s+from] + 
           0.5769*workHil[r-2][s+from] - 
           0.5769*workHil[r-4][s+from] - 
           0.0962*workHil[r-6][s+from]) * (0.075*workHil[r-1][s+_period] + 0.54));
}