Normalized, Variety, Fast Fourier Transform Explorer [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("Normalized, Variety, Fast Fourier Transform Explorer [Loxx]",
     shorttitle = "NVFFTE [Loxx]",
     overlay = false,
     max_lines_count = 500,
     precision = 8)

greencolor = #2DD204
redcolor = #D2042D

_fastsinetransform(float[] a, int tnn, bool inversefst)=>
    int j = 0
    int tm = 0
    int n2 = 0
    float sum = 0.
    float y1 = 0.
    float y2 = 0.
    float theta = 0.
    float wi = 0.
    float wr = 0.
    float wpi = 0.
    float wpr = 0.
    float wtemp = 0.
    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 istep = 0
    int isign = 0
    float tempr = 0.
    float tempi = 0.

    theta := math.pi/tnn
    wr := 1.0
    wi := 0.0
    wpr :=  - 2.0 * math.pow(math.sin(0.5 * theta),2)
    wpi := math.sin(theta)
    array.set(a, 0, 0.)
    tm := tnn / 2
    n2 := tnn + 2
    for jx = 2 to tm  +  1
        wtemp := wr
        wr := wr * wpr - wi * wpi + wr
        wi := wi * wpr + wtemp * wpi  +  wi
        y1 := wi * (array.get(a, jx - 1) + array.get(a, n2 - jx - 1))
        y2 := 0.5 * (array.get(a, jx - 1) - array.get(a, n2 - jx - 1))
        array.set(a, jx - 1, y1 + y2)
        array.set(a, n2 - jx - 1,  y1 - y2)
    ttheta := 2.0 * math.pi/tnn
    c1 := 0.5
    c2 := - 0.5
    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

    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))
    sum := 0.0
    array.set(a, 0, 0.5 * array.get(a, 0))
    array.set(a, 1, 0.)
    for jj = 0 to tm - 1
        j := 2 * jj + 1
        sum := sum + array.get(a, j - 1)
        array.set(a, j - 1, array.get(a, j))
        array.set(a, j, sum)
    if inversefst
        for jx = 1 to tnn
            array.set(a, jx - 1, array.get(a, jx - 1) * 2 / tnn)
    a

_fastcosinetransform(float[] a, int tnn, bool inversefct)=>
    int j = 0
    int n2 = 0
    float sum = 0.
    float y1 = 0.
    float y2 = 0.
    float theta = 0.
    float wi = 0.
    float wpi = 0.
    float wr = 0.
    float wpr = 0.
    float wtemp = 0.
    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 istep = 0
    int isign = 0
    float tempr = 0.
    float tempi = 0.

    while true
        if (tnn == 1)
            y1 := array.get(a, 0)
            y2 := array.get(a, 1)
            array.set(a, 0, 0.5 * (y1 + y2))
            array.set(a, 1, 0.5 * (y1 - y2))
            if inversefct
                array.set(a, 0, array.get(a, 0) * 2)
                array.set(a, 1, array.get(a, 1) * 2)
                break
        wi := 0
        wr := 1
        theta := math.pi / tnn
        wtemp := math.sin(theta * 0.5)
        wpr := -2.0 * wtemp * wtemp
        wpi := math.sin(theta)
        sum := 0.5 * (array.get(a, 0) - array.get(a, tnn))
        array.set(a, 0, 0.5 * (array.get(a, 0) + array.get(a, tnn)))
        n2 := tnn + 2
        for jx = 2 to tnn / 2
            wtemp := wr
            wr := wtemp * wpr - wi * wpi + wtemp
            wi := wi * wpr + wtemp * wpi + wi
            y1 := 0.5 * (array.get(a, jx - 1) + array.get(a, n2 - jx - 1))
            y2 := array.get(a, jx - 1) - array.get(a, n2 - jx - 1)
            array.set(a, jx - 1, y1 - wi * y2)
            array.set(a, n2 - jx - 1, y1 + wi * y2)
            sum := sum + wr * y2

        ttheta := 2.0 * math.pi / tnn
        c1 := 0.5
        c2 := -0.5
        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
        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))
        array.set(a, tnn, array.get(a, 1))
        array.set(a, 1, sum)
        j := 4
        while (j <= tnn)
            sum := sum + array.get(a, j - 1)
            array.set(a, j - 1, sum)
            j := j + 2
        if (inversefct)
            for jx = 0 to tnn
                array.set(a, jx, array.get(a, jx) * 2 / tnn)
        break
    a

_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

//Normalize data
_InSigNormalize(float[] aa)=>
    float sum_sqrt = 0.
    int element_count = array.size(aa)
    for i = 0 to element_count - 1
        sum_sqrt += math.pow(array.get(aa, i), 2)

    sum_sqrt := math.sqrt(sum_sqrt)
    if (sum_sqrt != 0)
        for i = 0 to element_count - 1
            array.set(aa, i, array.get(aa, i) / sum_sqrt)
    aa

src = input.source(open, "Source")
windowper = input.string("1024", "Window Period", options = ["16", "32", "64", "128", "256", "512", "1024", "2048"])
SS = input.int(20, "Smoothing Period")

type = input.string("Fast Sine", "Fast Fourier Transform Type", options =["Fast Cosine", "Fast Sine", "Real Fast"])

Fmin1 = input.int(3, "Minimum Frequency Filter")
Fmax1 = input.int(21, "Maximum Frequency Filter")
barsback = input.int(0, "Last Bar", minval = 0)
Inverse = input.bool(true, "Add inverse step?")

Window = str.tonumber(windowper)

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

var aa = array.new<float>(N + (type == "Fast Cosine" ? 1 : 0), 0.)

var pvlines = array.new_line(0)
countout = N >= 256 ? 250 : N

if barstate.isfirst
    for i = 0 to countout - 1
        array.push(pvlines, line.new(na, na, na, na))

if barstate.islast
    for i = 0 to N - 1
        array.set(aa, i, nz(src[i + barsback]))

    M = array.size(aa)
    end = M - 1
    ss = int(math.min(SS, M))

    _InSigNormalize(aa)

    if type == "Fast Cosine"
        _fastcosinetransform(aa, N, false)
    else if type == "Fast Sine"
        _fastsinetransform(aa, N, false)
    else
        _realFastFourierTransform(aa, N, false)

    for k = 0 to end
        if k >= ss
            array.set(aa, k, 0.)

    Fmax1 := math.min(Fmax1, array.size(aa) - 1)

    for i = 0 to array.size(aa) - 1
        if i < Fmin1 or i > Fmax1
            array.set(aa, i, 0.)

    if (Inverse)
        if type == "Fast Cosine"
            _fastcosinetransform(aa, N, true)
        else if type == "Fast Sine"
            _fastsinetransform(aa, N, true)
        else
            _realFastFourierTransform(aa, N, true)

    xm = array.copy(aa)

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

    while i < N  and i < array.size(xm) - 2

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

        i += skipper
        j += 1

plot(0, color = bar_index % 2 ? color.gray : na)