// This source code is subject to the terms of the Mozilla Public License 2.0 at https://mozilla.org/MPL/2.0/
// © yatrader2

// The Goertzel algorithm is a technique in digital signal processing (DSP) for efficient evaluation of the individual terms of the discrete Fourier transform (DFT).
// https://en.wikipedia.org/wiki/Goertzel_algorithm

// See also:
// Python implementation of the Goertzel algorithm for calculating DFT terms
// https://gist.github.com/sebpiq/4128537

// See also:
// https://www.st.com/resource/en/design_tip/dm00446805-the-goertzel-algorithm-to-compute-individual-terms-of-the-discrete-fourier-transform-dft-stmicroelectronics.pdf


//@version=4
study("FUNCTION: Goertzel algorithm / filter -- Calculate DFT of a specific frequency bin", overlay=false)

// CONSTANTS

pi=2*asin(1)
PI= pi
tau = 2*pi


// The following function takes a vector x (your source) made of N samples, and computes the k-th frequency coefficient (k from 0 to N-1).
// The real part is returned in as Inphase and the imaginary part is returned as Quadrature.
// The Magnitude and Theta are also calculated

// FUNCTIONS

// Goertzel algorithm for integer k
goertzel_int(x,k,N)=>
    // Because our signal is a complex signal we have to comput both the real and imaginary
    w = 2*pi*k/N     // omega = (2 pi * number of bars in period)/Window length
    w_real = cos(w) // cos of omega (the real portion of omega)
    w_imag = sin(w) // sin of omega (the complex portion of omega)
    c = 2*w_real    // Coefficient 2* cos of omega (real part)
    y  = 0.0
    z1 = 0.0
    z2 = 0.0
    for n = 0 to N-1
        y := x[n] + c*z1 - z2
        z2 := z1
        z1 := y
    Inphase    = w_real*z1 - z2   // Real
    Quadrature = w_imag*z1        // Imag
    Power  = pow(z2,2)+pow(z1,2)-  (c*z1*z2)
    Theta = atan(Inphase/Quadrature)*2*pi // In radians
    [Inphase, Quadrature, Power, Theta]


// To do implement hamming
hamm_f(_n,_N)=>
    hamm = 0.54 - 0.46 * cos(tau*(_n/(_N-1)))

// Inputs
src=input(close)

// recirpocal of frequency
len=input(20, title="Period length (recpirocal of frequency)", type=input.integer)
block_len=input(40, title="Measurment Block Length (N in goertzel notation)")
// dohamm=input(false)

k = block_len/len

[InPhase_real, Quadrature, Power, Theta]= goertzel_int(src,k,block_len)

plot (InPhase_real, title="InPhase_real", color=color.red, display=display.none)
plot (Quadrature, title="Quadrature", color=color.blue, display=display.none)
plot (Power, title="Power", color=color.green)
plot (Theta, title="Theta", color=color.orange, display=display.none)