Goertzel Cycle Period [Loxx]
// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © loxx

//@version=5
//@description Goertzel algorithm to extract cycle periods for use in adaptive indicators see:https://en.wikipedia.org/wiki/Goertzel_algorithm

indicator("Goertzel Cycle Period [Loxx]", overlay = false, shorttitle="GCP [Loxx]", timeframe="", timeframe_gaps=true, max_bars_back = 3000)

//@function Extract cycle period using Goertzel Algorithm
//@param src, float, source of cycle data to be extracted.
//@param per, int, maximum look back window
//@returns int, Goertzel Cycle pereiod 
_goertzel_cycle_period(float src, int per)=>

    float[] goeWork1 = array.new_float(per * 3, 0.) //Power array
    float[] goeWork2 = array.new_float(per * 3, 0.) //Theta array
    float[] goeWork3 = array.new_float(per * 3, 0.) //Store cylces
    float[] goeWork4 = array.new_float(per * 3, 0.) //Store ingest values
    
    float temp1 = nz(src[3 * per - 1]) //get first bar for calculation
    float temp2 = (src - temp1) / (3 * per - 1) // Real 
    float temp3 = 0. //Imaginary 
    
    for k = 3 * per - 1 to 1 
        float temp = src[k - 1] - (temp1 + temp2 * (3 * per - k))
        array.set(goeWork4, k, temp)
        array.set(goeWork3, k, 0.)
        
    for k = 2 to per 
        float w = 0., float x = 0., float y = 0.
        float z = math.pow(k, -1)
        
        temp1 := 2.0 * math.cos(2.0 * math.pi * z)
        
        //The first stage calculates an intermediate sequence, s[n]
        for i = 3 * per-1  to 1  
            w := temp1 * x - y + nz(array.get(goeWork4, i))
            y := x
            x := w

        // The second stage applies the following filter to y[n]
        // 1. The first filter stage can be observed to be a second-order IIR filter with a direct-form structure TBD
        // 2. The second-stage filter can be observed to be a FIR filter, since its calculations do not use any of its past outputs. TBD
        
        // calculate Real component 
        temp2 := x - y / 2.0 * temp1 
        if (temp2 == 0.0)
            temp2 := 0.0000001
            
        // calculate Imaginary component 
        temp3 := y * math.sin(2.0 * math.pi * z) 
        
        // calculate Power
        array.set(goeWork1, k, math.pow(temp2, 2) + math.pow(temp3, 2)) 
        
        // calculate Theta
        array.set(goeWork2, k, math.atan(temp3 / temp2))  
        if (temp2 < 0.0)
            array.set(goeWork2, k, nz(array.get(goeWork2, k)) + math.pi)
        else
            if (temp3 < 0.0) 
                array.set(goeWork2, k, nz(array.get(goeWork2, k)) + 2.0 * math.pi)
                
    // Parse cycles
    for k = 3 to per - 1
        if (nz(array.get(goeWork1, k)) > nz(array.get(goeWork1, k + 1)) and nz(array.get(goeWork1, k)) > nz(array.get(goeWork1, k-1)))
            array.set(goeWork3, k, k * math.pow(10, -4))
        else 
            array.set(goeWork3, k, 0.0)
    
    // count cycles
    int number_of_cycles = 0
    for i = 0 to per + 1  
        if (nz(array.get(goeWork3, i)) > 0.0) 
            number_of_cycles := number_of_cycles + 1
    number_of_cycles


src = input.source(close, "Source", group = "Basic Settings")
per = input.int(120, "Max Period", group = "Basic Settings")

cperout = _goertzel_cycle_period(src, per)
cperout := cperout < 1 ? 1 : cperout

plot(cperout, color = color.white, linewidth = 2)



















