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 functionMonday, 21 December 2015
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.
Subscribe to:
Post Comments (Atom)
No comments:
Post a Comment