John Ehler's Sinewave Indicator Code

A reader recently inquired about my use of this indicator and so below I provide my Octave C++ .oct  version that I have been using for the past few years.

DEFUN_DLD ( sinewave_indicator, args, nargout )
{
octave_value_list retval_list ;
int nargin = args.length () ;
int vec_length = args(0).length () ;

// check the input argument
if ( nargin != 1 )
   {
   error ("Invalid argument. Input is a single price vector.") ;
   return retval_list ;
   }

if ( vec_length < 50 )
   {
   error ("Invalid argument. Input is a single price vector.") ;
   return retval_list ;
   }

if ( error_state )
   {
   error ("Invalid argument. Input is a single price vector.") ;
   return retval_list ;
   }
// end of input checking

// inputs
ColumnVector price = args(0).column_vector_value () ;

// outputs
ColumnVector sinewave( vec_length ) ;
ColumnVector sinewave_lead_1( vec_length ) ;
ColumnVector smoothperiod_out( vec_length ) ;
ColumnVector dcphase_vec( vec_length ) ;
ColumnVector sumperiod( vec_length ) ;
ColumnVector sum_period( vec_length ) ;
ColumnVector deltaphase( vec_length ) ;

// Declarations for calculations of period, phase & sine wave measurements
ColumnVector smooth( vec_length ) ;
ColumnVector period( vec_length ) ;          
ColumnVector smoothperiod( vec_length ) ; 
ColumnVector detrender( vec_length ) ; 
ColumnVector Q1( vec_length ) ; 
ColumnVector I1( vec_length ) ; 
ColumnVector jI( vec_length ) ; 
ColumnVector jQ( vec_length ) ; 
ColumnVector I2( vec_length ) ;  
ColumnVector Q2( vec_length ) ;  
ColumnVector sI2( vec_length ) ; 
ColumnVector sQ2( vec_length ) ; 
ColumnVector Re( vec_length ) ;
ColumnVector Im( vec_length ) ; 
ColumnVector sRe( vec_length ) ; 
ColumnVector sIm( vec_length ) ;  
int dcperiod ; 
double realpart ;
double imagpart ;  
double dcphase ;
double sum_deltaphase ;
int count ;

// unrolled loop to fill the first 5 elements of above calculation vectors ( unrolled for speed optimisation )
sinewave(0) = 0.0 ; sinewave(1) = 0.0 ; sinewave(2) = 0.0 ; sinewave(3) = 0.0 ; sinewave(4) = 0.0 ;
sinewave_lead_1(0) = 0.0 ; sinewave_lead_1(1) = 0.0 ; sinewave_lead_1(2) = 0.0 ; sinewave_lead_1(3) = 0.0 ; sinewave_lead_1(4) = 0.0 ;
smoothperiod_out(0) = 0.0 ; smoothperiod_out(1) = 0.0 ; smoothperiod_out(2) = 0.0 ; smoothperiod_out(3) = 0.0 ; smoothperiod_out(4) = 0.0 ;
dcphase_vec(0) = 0.0 ; dcphase_vec(1) = 0.0 ; dcphase_vec(2) = 0.0 ; dcphase_vec(3) = 0.0 ; dcphase_vec(4) = 0.0 ;

smooth(0) = 0.0 ; smooth(1) = 0.0 ; smooth(2) = 0.0 ; smooth(3) = 0.0 ; smooth(4) = 0.0 ;
period(0) = 0.0 ; period(1) = 0.0 ; period(2) = 0.0 ; period(3) = 0.0 ; period(4) = 0.0 ;            
smoothperiod(0) = 0.0 ; smoothperiod(1) = 0.0 ; smoothperiod(2) = 0.0 ; smoothperiod(3) = 0.0 ; smoothperiod(4) = 0.0 ;
detrender(0) = 0.0 ; detrender(1) = 0.0 ; detrender(2) = 0.0 ; detrender(3) = 0.0 ; detrender(4) = 0.0 ; 
Q1(0) = 0.0 ; Q1(1) = 0.0 ; Q1(2) = 0.0 ; Q1(3) = 0.0 ; Q1(4) = 0.0 ;
I1(0) = 0.0 ; I1(1) = 0.0 ; I1(2) = 0.0 ; I1(3) = 0.0 ; I1(4) = 0.0 ; 
jI(0) = 0.0 ; jI(1) = 0.0 ; jI(2) = 0.0 ; jI(3) = 0.0 ; jI(4) = 0.0 ;
jQ(0) = 0.0 ; jQ(1) = 0.0 ; jQ(2) = 0.0 ; jQ(3) = 0.0 ; jQ(4) = 0.0 ;
I2(0) = 0.0 ; I2(1) = 0.0 ; I2(2) = 0.0 ; I2(3) = 0.0 ; I2(4) = 0.0 ; 
Q2(0) = 0.0 ; Q2(1) = 0.0 ; Q2(2) = 0.0 ; Q2(3) = 0.0 ; Q2(4) = 0.0 ;
sI2(0) = 0.0 ; sI2(1) = 0.0 ; sI2(2) = 0.0 ; sI2(3) = 0.0 ; sI2(4) = 0.0 ;
sQ2(0) = 0.0 ; sQ2(1) = 0.0 ; sQ2(2) = 0.0 ; sQ2(3) = 0.0 ; sQ2(4) = 0.0 ;
Re(0) = 0.0 ; Re(1) = 0.0 ; Re(2) = 0.0 ; Re(3) = 0.0 ; Re(4) = 0.0 ;
Im(0) = 0.0 ; Im(1) = 0.0 ; Im(2) = 0.0 ; Im(3) = 0.0 ; Im(4) = 0.0 ;
sRe(0) = 0.0 ; sRe(1) = 0.0 ; sRe(2) = 0.0 ; sRe(3) = 0.0 ; sRe(4) = 0.0 ;
sIm(0) = 0.0 ; sIm(1) = 0.0 ; sIm(2) = 0.0 ; sIm(3) = 0.0 ; sIm(4) = 0.0 ;
            
 for ( octave_idx_type ii (5) ; ii < vec_length ; ii++ ) // Start the main loop
     {
       
     // smooth the price for hilbert calculations
     smooth(ii) = (4.0 * price(ii) + 3.0 * price(ii-1) + 2.0 * price(ii-2) + price(ii-3) ) / 10.0 ; 
     
     // Detrend the input
     detrender(ii) = (0.0962 * smooth(ii) + 0.5769 * smooth(ii-2) - 0.5769 * smooth(ii-4) - 0.0962 * smooth(ii-6)) * (0.075 * period(ii-1) + 0.54) ;

     // Compute  InPhase and Quadrature components 
     Q1(ii) = (0.0962 * detrender(ii) + 0.5769 * detrender(ii-2) - 0.5769 * detrender(ii-4) - 0.0962 * detrender(ii-6)) * (0.075 * period(ii-1) + 0.54) ;
     I1(ii) = detrender(ii-3) ;

     // Advance the phase of  I1 and Q1 by 90 degrees
     jI(ii) = (0.0962 * I1(ii) + 0.5769 * I1(ii-2) - 0.5769 * I1(ii-4) - 0.0962 * I1(ii-6)) * (0.075 * period(ii-1) + 0.54) ;
     jQ(ii) = (0.0962 * Q1(ii) + 0.5769 * Q1(ii-2) - 0.5769 * Q1(ii-4) - 0.0962 * Q1(ii-6)) * (0.075 * period(ii-1) + 0.54) ;

     // Phasor addition for 3 bar averaging
     I2(ii) = I1(ii) - jQ(ii) ;
     Q2(ii) = Q1(ii) + jI(ii) ;

     // Smooth the  I and Q components before applying the discriminator
     sI2(ii) = 0.2 * I2(ii) + 0.8 * sI2(ii-1) ;
     sQ2(ii) = 0.2 * Q2(ii) + 0.8 * sQ2(ii-1) ;

     // Homodyne Discriminator
     Re(ii) = sI2(ii) * sI2(ii-1) + sQ2(ii) * sQ2(ii-1) ;
     Im(ii) = sI2(ii) * sQ2(ii-1) - sQ2(ii) * sI2(ii-1) ;
     sRe(ii) = 0.2 * Re(ii) + 0.8 * sRe(ii-1) ;
     sIm(ii) = 0.2 * Im(ii) + 0.8 * sIm(ii-1) ; 

       if ( (sIm(ii) > 0.0 || sIm(ii) < 0.0) && (sRe(ii) > 0.0 || sRe(ii) < 0.0) )
       { 
       period(ii) = 360.0 / ( ((atan(sIm(ii) / sRe(ii))) * 180.0) / PI ) ;
       }
       else
       {
       period(ii) = period(ii-1) ;
       }

       if ( period(ii) > 1.5 * period(ii-1) )
       {
       period(ii) = 1.5 * period(ii-1) ;
       }

       if ( period(ii) < 0.67 * period(ii-1) )
       {
       period(ii) = 0.67 * period(ii-1) ;
       }

       if ( period(ii) < 6.0 )
       {
       period(ii) = 6.0 ;
       }

       if ( period(ii) > 50.0 )
       {
       period(ii) = 50.0 ;
       }
 
     period(ii) = 0.2 * period(ii) + 0.8 * period(ii-1) ;
     smoothperiod(ii) = 0.33 * period(ii) + 0.67 * smoothperiod(ii-1) ;

     // Compute Dominant Cycle
     dcperiod = int ( smoothperiod(ii) + 0.5 ) ;
     realpart = 0.0 ;
     imagpart = 0.0 ;
     dcphase = 0.0 ;

      for ( octave_idx_type jj (0) ; jj <= ( dcperiod - 1 ) ; jj++ )
          {
          realpart += sin( PI/180.0 * 360.0 * jj / dcperiod ) * ( smooth(ii-jj) ) ;
          imagpart += cos( PI/180.0 * 360.0 * jj / dcperiod ) * ( smooth(ii-jj) ) ;
          }

      if ( fabs( imagpart ) > 0.0 )
         {
         dcphase = atan( realpart / imagpart ) * 180.0 / PI ;
         }
      else if ( fabs( imagpart ) < 0.001 )
              {
              if ( realpart < 0.0 )
                 {
                 dcphase -= 90.0 ;
                 }
              else if ( realpart > 0.0 )
                      {
                      dcphase += 90.0 ;
                      }
              }
         dcphase += 90.0 ;

      // Compensate for one bar lag of the 4 bar weighted moving average
      dcphase += 360.0 / smoothperiod(ii) ;

      if ( imagpart < 0.0 )
         dcphase += 180.0 ;

      if ( dcphase > 315.0 )
         dcphase -= 360.0 ;
     
     // phase output 
     dcphase_vec(ii) = dcphase ;
     
     //Now compute a differential phase, resolve phase wraparound, and limit delta phase errors
     deltaphase(ii) = dcphase_vec(ii) - dcphase_vec(ii-1) ;
     
     if ( dcphase_vec(ii-1) > 270.0 && dcphase_vec(ii) < 90.0 )
        {
 deltaphase(ii) = 360.0 - dcphase_vec(ii-1) + dcphase_vec(ii) ;
 }
 
     if ( deltaphase(ii) < 1.0 )
        { 
 deltaphase(ii) = 1.0 ;
 }
 
     if ( deltaphase(ii) > 60.0 )
        {
 deltaphase(ii) = 60.0 ;
 }

     // Sum Deltaphases to reach 360 degrees. The sum is the instantaneous period.
     sum_period(ii) = 0.0 ;
     sum_deltaphase = 0.0 ;
     count = 0 ;
     
     while ( sum_deltaphase < 360.0 ) 
           {
           sum_deltaphase += deltaphase(ii-count) ;
           count ++ ;
    sum_period(ii) = count ;
           } 

      // Resolve Instantaneous Period errors and smooth
      if ( sum_period(ii) == 0.0 )
         {
  sum_period(ii) = sum_period(ii-1) ;
  }
  
      sumperiod(ii) = 0.25 * sum_period(ii) + 0.75 * sum_period(ii-1) ;

     // sinewave output
     sinewave(ii) = sin( dcphase * PI / 180.0 ) ;
     
     // one bar leading function
     sinewave_lead_1(ii) = sin( ( dcphase + 360.0 / smoothperiod(ii) ) * PI / 180.0 ) ;
     
     // period output
     smoothperiod_out(ii) = floor ( smoothperiod(ii) + 0.5 ) ;
        
     } // end of main ii loop
      
 retval_list(3) = dcphase_vec ;
 retval_list(2) = smoothperiod_out ;
 retval_list(1) = sinewave_lead_1 ;                                                                  
 retval_list(0) = sinewave ;

return retval_list ; 
                                                                       
} // end of function

This is a straightforward conversion of the code available from here. A nice intro to how it can be used is here and Ehler's own website can be found here. 

Denoising Autoencoder MATLAB/Octave Code

Following on from my last post I have been looking for Octave code for the denoising autoencoder to avoid reinventing the wheel and writing it myself from scratch, and luckily I have found two options. The first is a tutorial on autoencoders, by a Piotr Mirowski, which has a link to a Github page with code. The second is a toolbox for dimensionality reduction, by Laurens van der Maaten, which has autoencoders as one of its options.

As of now I'm not sure which one I'll use, and perhaps I might yet write my own code heavily reusing elements from both of the above. More in due course.