Real-Fast Fourier Transform of Price w/ Linear Regression [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
indicator("Real-Fast Fourier Transform of Price w/ Linear Regression [Loxx]", 
     shorttitle = "RFFTP [Loxx]",
     overlay = true, 
     max_lines_count = 500)
     
greencolor = #2DD204
redcolor = #D2042D 

//+-------------------------------------------------------------------------+
//| Real Fast Fourier Transform |
///+-------------------------------------------------------------------------+
//The algorithm performs a fast Fourier transform of a real
//function defined by n samples on the real axis.
//
//Depending on the passed parameters, it can be executed
//both direct and inverse conversion.
//
// Input parameters:
//    tnn - Number of function values. Must be degree
//            deuces. Algorithm does not validate
//            passed value.
//    a - array [0 .. nn-1] of Real
//            Function values.
//    InverseFFT
//        - the direction of the transformation.
//            True if reverse, False if direct.
//         
// Output parameters:
//    a - the result of the transformation. 
//
// For more details, see description on the site:
// https://www.alglib.net/fasttransforms/fft.php 
// https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.205.4523&rep=rep1&type=pdf
//+-------------------------------------------------------------------------+

_realFastFourierTransform(float[] a, int tnn, bool inversefft)=>
    float twr = 0.
    float twi = 0.
    float twpr = 0.
    float twpi = 0.
    float twtemp = 0.
    float ttheta = 0.
    int i = 0
    int i1 = 0
    int i2 = 0
    int i3 = 0
    int i4 = 0
    float c1 = 0.
    float c2 = 0.
    float h1r = 0.
    float h1i = 0.
    float h2r = 0.
    float h2i = 0.
    float wrs = 0.
    float wis = 0.
    int nn = 0
    int n = 0
    int mmax = 0
    int m = 0
    int j = 0
    int istep = 0
    int isign = 0
    float wtemp = 0.
    float wr = 0.
    float wpr = 0.
    float wpi = 0.
    float wi = 0.
    float theta = 0.
    float tempr = 0.
    float tempi = 0.

    if (tnn != 1) 
        if (not inversefft) 
            ttheta := 2.0 * math.pi / tnn
            c1 := 0.5
            c2 := -0.5
        else 
            ttheta := 2.0 * math.pi / tnn
            c1 := 0.5
            c2 := 0.5
            ttheta := -ttheta
            twpr := -2.0 * math.pow(math.sin(0.5 * ttheta), 2)
            twpi := math.sin(ttheta)
            twr := 1.0 + twpr
            twi := twpi
            
            for ix = 2 to tnn / 4 + 1 
                i1 := ix + ix - 2
                i2 := i1 + 1
                i3 := tnn + 1 - i2
                i4 := i3 + 1
                wrs := twr
                wis := twi
                h1r := c1 * (array.get(a, i1) + array.get(a, i3))
                h1i := c1 * (array.get(a, i2) - array.get(a, i4))
                h2r := -c2 * (array.get(a, i2) + array.get(a, i4))
                h2i := c2 * (array.get(a, i1) - array.get(a, i3))
                array.set(a, i1,  h1r + wrs * h2r - wis * h2i)
                array.set(a, i2,  h1i + wrs * h2i + wis * h2r)
                array.set(a, i3,  h1r - wrs * h2r + wis * h2i)
                array.set(a, i4, -h1i + wrs * h2i + wis * h2r)
                twtemp := twr
                twr := twr * twpr - twi * twpi + twr
                twi := twi * twpr + twtemp * twpi + twi
            
            
            h1r := array.get(a, 0)
            array.set(a, 0, c1 * (h1r + array.get(a, 1)))
            array.set(a, 1, c1 * (h1r - array.get(a, 1)))
        
        if (inversefft) 
            isign := -1
        else 
            isign := 1
        
        n := tnn
        nn := tnn / 2
        j := 1
        for ii = 1 to nn  
            i := 2 * ii - 1
            if (j > i) 
                tempr := array.get(a, j - 1)
                tempi := array.get(a, j)
                array.set(a, j - 1, array.get(a, i - 1))
                array.set(a, j, array.get(a, i))
                array.set(a, i - 1, tempr)
                array.set(a, i, tempi)
            
            m := n / 2
            while (m >= 2 and j > m) 
                j := j - m
                m := m / 2
            
            j := j + m
        
        mmax := 2
        while (n > mmax) 
            istep := 2 * mmax
            theta := 2.0 * math.pi / (isign * mmax)
            wpr := -2.0 * math.pow(math.sin(0.5 * theta), 2)
            wpi := math.sin(theta)
            wr := 1.0
            wi := 0.0
            for ii = 1 to mmax / 2  
                m := 2 * ii - 1
                for jj = 0 to (n - m) / istep 
                    i := m + jj * istep
                    j := i + mmax
                    tempr := wr * array.get(a, j - 1) - wi * array.get(a, j)
                    tempi := wr * array.get(a, j) + wi * array.get(a, j - 1)
                    array.set(a, j - 1, array.get(a, i - 1) - tempr)
                    array.set(a, j, array.get(a, i) - tempi)
                    array.set(a, i - 1, array.get(a, i - 1) + tempr)
                    array.set(a, i, array.get(a, i) + tempi)
                
                wtemp := wr
                wr := wr * wpr - wi * wpi + wr
                wi := wi * wpr + wtemp * wpi + wi
            
            mmax := istep
    
        if (inversefft) 
            for ix = 1 to 2 * nn
                array.set(a, ix - 1, array.get(a, ix - 1) / nn)

        if (not inversefft) 
            twpr := -2.0 * math.pow(math.sin(0.5 * ttheta), 2)
            twpi := math.sin(ttheta)
            twr := 1.0 + twpr
            twi := twpi
            
            for ix = 2 to tnn / 4 + 1 
                i1 := ix + ix - 2
                i2 := i1 + 1
                i3 := tnn + 1 - i2
                i4 := i3 + 1
                wrs := twr
                wis := twi
                h1r := c1 * (array.get(a, i1) + array.get(a, i3))
                h1i := c1 * (array.get(a, i2) - array.get(a, i4))
                h2r := -c2 * (array.get(a, i2) + array.get(a, i4))
                h2i := c2 * (array.get(a, i1) - array.get(a, i3))
                array.set(a, i1, h1r + wrs * h2r - wis * h2i)
                array.set(a, i2, h1i + wrs * h2i + wis * h2r)
                array.set(a, i3, h1r - wrs * h2r + wis * h2i)
                array.set(a, i4, -h1i + wrs * h2i + wis * h2r)
                twtemp := twr
                twr := twr * twpr - twi * twpi + twr
                twi := twi * twpr + twtemp * twpi + twi

            h1r := array.get(a, 0)
            array.set(a, 0, h1r + array.get(a, 1))
            array.set(a, 1, h1r - array.get(a, 1))
    a

//linear regression value
_iLRValue(float[] src, float period, float slope)=>
    float SumXY = 0.
    float SumY = 0.
    float sumX = period * (period - 1) / 2.0
    float sumXX = period * (period - 1) * (2.0 * period - 1.0) / 6.0
    float divisor = math.pow(sumX, 2.0) - period * sumXX

    slope1 = slope
    for k = 0 to period - 1
        SumXY += k * array.get(src, k)
        SumY += array.get(src, k)
  
    if (divisor != 0)
        slope1 := (period * SumXY - sumX * SumY) / divisor
    else 
        slope1 := 0
    out = ((SumY - slope1 * sumX) / period + slope1 * (period - 1.0))
    [out, slope1]

src = input.source(close, "Source", group = "Basic Settings")
uselreg = input.bool(true, "Use Linear Regression?", group = "Basic Settings")
Windowin = input.string("512", "Window", options = ["4", "8", "16", "32", "64", "128", "256", "512", "1024", "2048"], group = "Basic Settings")
Treshold = input.float(0.03, "Treshold", step = 0.001, group = "Basic Settings")
dtrendper = input.int(0, "Detrended Line Period", group = "Basic Settings")
barsback = input.int(30, "Last Bar", minval = 0, group = "Basic Settings")
mutebars = input.bool(true, "Bar color muting?", group = "Basic Settings")

//+------------------------------------------------------------------+
//| Initialize variables
//+------------------------------------------------------------------+

var pvlines = array.new_line(0)
var lrlines = array.new_line(0)

//+------------------------------------------------------------------+
//| Force input Window period to lowest nearest power of 2,
//| we alreadsy solved to make sure user can 
//| only select powers of 2
//+------------------------------------------------------------------+

Window = str.tonumber(Windowin) 
n = int(math.log(Window) / math.log(2))
N = int(math.max(math.pow(2, n), 16))

barcolor(mutebars ? last_bar_index - bar_index >= barsback and  last_bar_index - bar_index  < N  ?  color.gray : na : na)

//+------------------------------------------------------------------+
//| Initialize lines for drawing later
//+------------------------------------------------------------------+
countout = Window >= 256 ? 250 : Window
if barstate.isfirst
    for i = 0 to countout - 1
        array.push(pvlines, line.new(na, na, na, na))
        array.push(lrlines, line.new(na, na, na, na))

//+------------------------------------------------------------------+
//| Core compute 
//+------------------------------------------------------------------+
if barstate.islast

    //+------------------------------------------------------------------+
    //| Initial variables  
    //+------------------------------------------------------------------+
    dlLength = dtrendper > N or dtrendper < 1 ? N :  dtrendper
    work = array.new<float>(N, 0.)
    price = array.new<float>(N, 0.)
    lr = array.new<float>(N, 0.)
    fft = array.new<float>(N, 0.)
    
    for i = 0 to N - 1
        array.set(price, i, nz(src[i + barsback]))
        
    endValue = 0.
    slope = 0.
    
    //+------------------------------------------------------------------+
    //| Apply Linear Regression and detrending
    //+------------------------------------------------------------------+
    if (uselreg)
        [dv, slp] = _iLRValue(price, dlLength, slope)
        endValue := dv
        slope := slp
    else 
        endValue := array.get(price, 0)
        slope := (endValue - array.get(price, dlLength - 1)) / dlLength
    
    for i = N - 1 to 0 
        array.set(lr, i, endValue - slope * i)
        array.set(work, i, array.get(price, i) - array.get(lr, i))
    
    //+------------------------------------------------------------------+
    //| Apply Regular FFT and calculate powers and modified by threshhold
    //+------------------------------------------------------------------+
    _realFastFourierTransform(work, N, false)
    
    for i = 1 to N / 2 - 1
        amp = math.sqrt(math.pow(array.get(work, i * 2), 2) + math.pow(array.get(work, i * 2 + 1), 2))
        add = 0.
        if (not (array.get(work, i * 2) > 0)) 
            add := math.pi
        phase = math.atan(array.get(work, i * 2 + 1) / array.get(work, i * 2))
        array.set(work, i * 2, amp * math.cos(phase + add))
        array.set(work, i * 2 + 1, amp * math.sin(phase + add))
        
    power = 0.
    for i = 1 to N / 2 - 1
        power += math.sqrt(math.pow(array.get(work, i * 2), 2) + math.pow(array.get(work, i * 2 + 1), 2))
    treshold = power * Treshold
    
    for i = 1 to  N / 2 - 1
        tmp = math.sqrt(math.pow(array.get(work, i * 2), 2) + math.pow(array.get(work, i * 2 + 1), 2))
        if (tmp < treshold) 
            array.set(work, i * 2, 0)
            array.set(work, i * 2 + 1, 0)
    
    //+------------------------------------------------------------------+
    //| Apply Inverse FFT and final output
    //+------------------------------------------------------------------+
    _realFastFourierTransform(work, N , true)
    for i = 0 to N - 1
        array.set(fft, i, array.get(work, i) + array.get(lr, i))

    skipper = Window == 2048 ? 8 : Window == 1024 ? 4 : Window == 512 ? 2 : 1
    int i = 0
    int j = 0

    //+------------------------------------------------------------------+
    //| Draw lines w/ skipping to stay within 500 line limit
    //+------------------------------------------------------------------+
    while i < N and i < array.size(fft) - 2 

        if j > array.size(pvlines) - 1
            break
        pvline = array.get(pvlines, j)
        colorout = i < array.size(fft) - 2 ? array.get(fft, i) > array.get(fft, i + skipper) ? greencolor : redcolor : na
        line.set_xy1(pvline, bar_index - i - skipper - barsback, array.get(fft, i + skipper))
        line.set_xy2(pvline, bar_index - i - barsback, array.get(fft, i))
        line.set_color(pvline, colorout)
        line.set_style(pvline, line.style_solid)
        line.set_width(pvline, 5)

        lrline = array.get(lrlines, j)
        colorout2 = i < array.size(lr) - 2 ? array.get(lr, i) > array.get(lr, i + skipper) ? greencolor : redcolor : na
        line.set_xy1(lrline, bar_index - i - skipper - barsback, array.get(lr, i + skipper))
        line.set_xy2(lrline, bar_index - i - barsback, array.get(lr, i))
        line.set_color(lrline, colorout2)
        line.set_style(lrline, line.style_solid)
        line.set_width(lrline, 2)
        
        i += skipper
        j += 1

