Giving Up on Recursive Sine Formula for Period Calculation

I have spent the last few weeks trying to get my recursive sine wave formula for period calculations to work, but try as I might I can only get it to do so under ideal theoretical conditions. Once any significant noise, trend or combination thereof is introduced the calculations explode and give meaningless results. In light of this, I am no longer going to continue this work.

Apart from the above work I have also been doing my usual online research and have come across John Ehler's autocorrelation periodogram for period measurement, and below is my Octave C++ .oct implementation of it.

DEFUN_DLD ( autocorrelation_periodogram, args, nargout,
"-*- texinfo -*-\n\
@deftypefn {Function File} {} autocorrelation_periodogram (@var{input_vector})\n\
This function takes an input vector ( price ) and outputs the dominant cycle period,\n\
calculated from the autocorrelation periodogram spectrum.\n\
@end deftypefn" )

{
octave_value_list retval_list ;
int nargin = args.length () ;

// check the input arguments
if ( nargin != 1 ) // there must be a price vector only
   {
   error ("Invalid arguments. Input is a price vector only.") ;
   return retval_list ;
   }

if ( args(0).length () < 4 )
   {
   error ("Invalid argument length. Input is a price vector of length >= 4.") ;
   return retval_list ;
   }

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

ColumnVector input = args(0).column_vector_value () ;
ColumnVector hp = args(0).column_vector_value () ; hp.fill( 0.0 ) ;
ColumnVector smooth = args(0).column_vector_value () ; smooth.fill( 0.0 ) ;
ColumnVector corr ( 49 ) ; corr.fill( 0.0 ) ;
ColumnVector cosine_part ( 49 ) ; cosine_part.fill( 0.0 ) ;
ColumnVector sine_part ( 49 ) ; sine_part.fill( 0.0 ) ;
ColumnVector sq_sum ( 49 ) ; sq_sum.fill( 0.0 ) ;
ColumnVector R1 ( 49 ) ; R1.fill( 0.0 ) ;
ColumnVector R2 ( 49 ) ; R2.fill( 0.0 ) ;  
ColumnVector pwr ( 49 ) ; pwr.fill( 0.0 ) ;
ColumnVector dominant_cycle = args(0).column_vector_value () ; dominant_cycle.fill( 0.0 ) ;  
   
double avglength = 3.0 ;
double M ;
double X ; double Y ;   
double Sx ; double Sy ; double Sxx ; double Syy ; double Sxy ;
double denom ;
double max_pwr = 0.0 ;
double Spx ; double Sp ;  

// variables for highpass filter, hard coded for a high cutoff period of 48 bars and low cutoff of 10 bars
double high_cutoff = 48.0 ; double low_cutoff = 10.0 ;  
double alpha_1 = ( cos( 0.707 * 2.0 * PI / high_cutoff ) + sin( 0.707 * 2.0 * PI / high_cutoff ) - 1.0 ) / cos( 0.707 * 2.0 * PI / high_cutoff ) ;   
double beta_1 = ( 1.0 - alpha_1 / 2.0 ) * ( 1.0 - alpha_1 / 2.0 ) ;
double beta_2 = 2.0 * ( 1.0 - alpha_1 ) ;
double beta_3 = ( 1.0 - alpha_1 ) * ( 1.0 - alpha_1 ) ;
 
// variables for super smoother
double a1 = exp( -1.414 * PI / low_cutoff ) ;
double b1 = 2.0 * a1 * cos( 1.414 * PI / low_cutoff ) ;
double c2 = b1 ;
double c3 = -a1 * a1 ;
double c1 = 1.0 - c2 - c3 ;
  
// calculate the automatic gain control factor, K
double K = 0.0 ;
double accSlope = -1.5 ; //acceptableSlope = 1.5 dB
double halfLC = low_cutoff / 2.0 ;
double halfHC = high_cutoff / 2.0 ;
double ratio = pow( 10 , accSlope / 20.0 ) ;
  
 if( halfHC - halfLC > 0.0 )
    {
  K = pow( ratio , 1.0 / ( halfHC - halfLC ) ) ;
    }

// loop to initialise hp and smooth
for ( octave_idx_type ii ( 2 ) ; ii < 49 ; ii++ ) // main loop
    {  
    // highpass filter components whose periods are < 48 bars
    hp(ii) = beta_1 * ( input(ii) - 2.0 * input(ii-1) + input(ii-2) ) + beta_2 * hp(ii-1) - beta_3 * hp(ii-2) ;
    
    // smooth with a super smoother filter
    smooth(ii) = c1 * ( hp(ii) + hp(ii-1) ) / 2.0 + c2 * smooth(ii-1) + c3 * smooth(ii-2) ;
    } // end of initial loop

for ( octave_idx_type ii ( 49 ) ; ii < args(0).length () ; ii++ ) // main loop
    {  
    // highpass filter components whose periods are < 48 bars
    hp(ii) = beta_1 * ( input(ii) - 2.0 * input(ii-1) + input(ii-2) ) + beta_2 * hp(ii-1) - beta_3 * hp(ii-2) ;
    
    // smooth with a super smoother filter
    smooth(ii) = c1 * ( hp(ii) + hp(ii-1) ) / 2.0 + c2 * smooth(ii-1) + c3 * smooth(ii-2) ;
      
      // Pearson correlation for each value of lag
      for ( octave_idx_type lag (0) ; lag <= high_cutoff ; lag++ ) 
          {
          // set the averaging length as M
          M = avglength ;
          if ( avglength == 0) 
             {
              M = double( lag ) ; 
             }
             
          Sx = 0.0 ; Sy = 0.0 ; Sxx = 0.0 ; Syy = 0.0 ; Sxy = 0.0 ;
            
            for ( octave_idx_type count (0) ; count < M - 1 ; count++ )
                {
                 X = smooth(ii-count) ; Y = smooth(ii-(lag+count)) ; 
                 Sx += X ; 
                 Sy += Y ;
                 Sxx += X * X ;
                 Sxy += X * Y ;
                 Syy += Y * Y ;
                }
             
            denom = ( M * Sxx - Sx * Sx ) * ( M * Syy - Sy * Sy ) ;    
            if ( denom > 0.0 )
               {
                corr(lag) = ( M * Sxy - Sx * Sy ) / sqrt( denom ) ;  
               }    
            
          } // end of Pearson correlation loop        
/*
    The DFT is accomplished by correlating the autocorrelation at each value of lag with the cosine and sine of each period of interest. 
    The sum of the squares of each of these values represents the relative power at each period.
*/                  
      for ( octave_idx_type period (low_cutoff) ; period <= high_cutoff ; period++ )
          {
           cosine_part( period ) = 0.0 ; sine_part( period ) = 0.0 ;
            
            for ( octave_idx_type N (3) ; N <= high_cutoff ; N++ )
                {
                 cosine_part( period ) += corr( N ) * cos( 2.0 * PI * double( N ) / double( period ) ) ;  
                 sine_part( period ) += corr( N ) * sin( 2.0 * PI * double( N ) / double( period ) ) ; 
                } // end of N loop
                
            sq_sum( period ) = cosine_part( period ) * cosine_part( period ) + sine_part( period ) * sine_part( period ) ;
                
          } // end of first period loop 
 
      // EMA is used to smooth the power measurement at each period          
      for ( octave_idx_type period (low_cutoff) ; period <= high_cutoff ; period++ )
          {
           R2( period ) = R1( period ) ;
           R1( period ) = 0.2 * sq_sum( period ) * sq_sum( period ) + 0.8 * R2( period ) ; 
          } // end of second period loop
      
    // Find maximum power level for normalisation      
    max_pwr = 0.0 ;      
          
      for ( octave_idx_type period (low_cutoff) ; period <= high_cutoff ; period++ )
          {
           if ( R1( period ) > max_pwr )
              {
                max_pwr = K * R1( period ) ;
              }
          } // end of third period loop
    
    // normalisation of power      
      for ( octave_idx_type period (low_cutoff) ; period <= high_cutoff ; period++ )
          {
           pwr( period ) = R1( period ) / max_pwr ;
          } // end of fourth period loop 
      
    // compute the dominant cycle using the centre of gravity of the spectrum
    Spx = 0.0 ; Sp = 0.0 ;
    
      for ( octave_idx_type period (low_cutoff) ; period <= high_cutoff ; period++ )
          {
           if ( pwr( period ) >= 0.5 )
              {
               Spx += double( period ) * pwr( period ) ;
               Sp += pwr( period ) ;  
              }
          } // end of fifth period loop
    
    if ( Sp != 0.0 )
       {
        dominant_cycle(ii) = Spx / Sp ;
       }      
          
    } // end of main loop
    
retval_list( 0 ) = dominant_cycle ;

return retval_list ;

} // end of function 

When applied directly to a theoretical but noisy sine wave series with a trend I find that this autocorrelation method performs better than my current period measurement algo, but on detrended data it is not as good. Since it is trivial to detrend price data, for now I am going to stick with my current method.