//------------------------------------------------------------------------------
//
//  Formula Name:    Ehler's filters and indicators
//  Author/Uploader: elvf 
//  E-mail:          
//  Date/Time Added: 2005-08-07 20:35:13
//  Origin:          US
//  Keywords:        Ehler's filters and indicators
//  Level:           medium
//  Flags:           function
//  Formula URL:     http://www.amibroker.com/library/formula.php?id=541
//  Details URL:     http://www.amibroker.com/library/detail.php?id=541
//
//------------------------------------------------------------------------------
//
//  I put together several functions, implementing Ehler's filters
//
//  and indicators into one include file.
//
//------------------------------------------------------------------------------

/*
All of the formulas below are from John Ehlers
*/
function ellipticFilter(P)
{
s[0]=P[0];
s[1]=P[1];
s[2]=P[2];
s[3]=P[3];

for(i=4;i<BarCount;i++)
{
s[i] = 0.13785*P[i] + 0.0007*P[i-1] + 0.13785*P[i-2] + 
       1.2103*s[i-1] - 0.4867*s[i-2];
}
return s;
}

function ellipticFilterZL(P)
{
// first, we need estimate of the velocity
Velocity=((P-Ref(P,-1))*2+Ref((P-Ref(P,-1)),-1))/3;

s[0]=P[0];
s[1]=P[1];
s[2]=P[2];
s[3]=P[3];

for(i=4;i<BarCount;i++)
{
s[i] = 0.13785*(P[i]+Velocity[i]) + 
       0.0007*(P[i-1]+Velocity[i-1]) + 
       0.13785*(P[i-2]+Velocity[i-2]) + 
       1.2103*s[i-1] - 0.4867*s[i-2];
}
return s;
}

function detrend(P)
{
d=P;
d2=P*0;
dt=P*0;
for(i=13;i<BarCount;i++)
{
 d[i] = 0.912*P[i]+0.088*d[i-6];
 d2[i] = (d[i] - d[i-6]) + (1.2*d2[i-6] - 0.7*d2[i-12]); 
 dt[i]=(d2[i-12]-d2[i-6])+(d2[i]-d2[i-6]);
}

return dt;
}

function ampl(P,range)
{
Imult=0.635;
Qmult=0.338;
InPhase[0]=0;
Quadrature[0]=0;
Amplitude[0]=0;

Value1=P*0;
Value2=Value1;

Inphase=P*0;
Quadrature=Inphase+1;

//Detrend Price
for(i=7;i<BarCount;i++)
{
Value1[i] = P[i] - P[i-7];

//Compute Hilbert Transform outputs
Inphase[i] = 1.25*(Value1[i-4] - Imult*Value1[i-2]) + Imult*InPhase[i-3];
Quadrature[i] = Value1[i-2] - Qmult*Value1[i] + Qmult*Quadrature[i-2];

//Compute smoothed Signal amplitude
Value2[i] = 0.2*(InPhase[i]*InPhase[i] + Quadrature[i]*Quadrature[i]) + 0.8*Value2[i-1];
}

//Compute smoothed SNR in Decibels, guarding against a divide by zero error, AND compensating for Filter loss
Value2=Max(Value2 ,0.001);
Amplitude = 0.25*(10*log(Value2/(Range*Range+0.0001))/log(10) + 1.9) + 0.75*Ref(Amplitude,-1);

return Amplitude;
}

// Ehlers Dominant Cycle Period
// from Ehlers, John F. Cybernetic Analysis for Stocks and Futures. Wiley. 2004. 
// Chapter 9, p. 107. Code on p. 111.

function CyclePeriod(array, alpha)
// Figure 9.4 on p. 111
{
  smooth = (array + 2*Ref(array, -1) + 2*Ref(array, -2) + Ref(array, -3))/6;

//  for(i = 0; i < 7; i++) cycle[i]=array[i]; // Initialize early values and as array

  for(i = 0; i < 6; i++) 
  {
     InstPeriod[i] = 0; // Initialize early values and as array
     DeltaPhase[i] = 0;
     cycle[i]=0;
     Period[i]=0;
  }

  for(i = 6; i < BarCount; i++)
  {
     cycle[i] = (1 - .5*alpha)*(1 - .5*alpha)*(smooth[i] - 2*smooth[i-1] + smooth[i-2]) +
                2*(1 - alpha)*cycle[i-1] - (1 - alpha)*(1 - alpha)*cycle[i-2];
     Q1[i] = (.0962*cycle[i] + .5769*cycle[i-2] -.5769*cycle[i-4] - .0962*cycle[i-6])*(.5 + .08*InstPeriod[i-1]);
     I1[i] = cycle[i-3];

     if(Q1[i] != 0 AND Q1[i-1] != 0) 
        DeltaPhase[i] = (I1[i]/Q1[i] - I1[i-1]/Q1[i-1])/(1 + I1[i]*I1[i-1]/(Q1[i]*Q1[i-1]));
     if(DeltaPhase[i] < 0.1) DeltaPhase[i] = 0.1;
     if(DeltaPhase[i] > 1.1) DeltaPhase[i] = 1.1;

     //----- Speed up the median calculation by placing it inline

     mlen = 5;
     for(k = mlen - 1; k >= 0; k--) {temparray[k] = DeltaPhase[i + k - (mlen - 1)];}

     temp=0;
     for(k = mlen - 1; k > 0; k--)
     {for (j = mlen - 1; j > 0; j--)
       {if (temparray[j-1] > temparray[j])
         {
           temp = temparray[j-1];
           temparray[j-1] = temparray[j];
           temparray[j] = temp;
         }
       }
     }
     MedianDelta[i] = temparray[mlen - 1 - (mlen / 2)];

     //----- End median calculation

     if(MedianDelta[i] == 0) DC[i] = 15; 
     else DC[i] = 6.28318/MedianDelta[i] + .5;

     InstPeriod[i] = .33*DC[i] + .67*InstPeriod[i-1];
     Period[i] = .15*InstPeriod[i] + .85*Period[i-1];
  }
  return Period;
}