loxxfft
// 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 This code is a library for performing Fast Fourier Transform (FFT) operations. FFT is an algorithm that can quickly compute the discrete Fourier transform (DFT) of a sequence. The library includes functions for performing FFTs on both real and complex data. It also includes functions for fast correlation and convolution, which are operations that can be performed efficiently using FFTs. Additionally, the library includes functions for fast sine and cosine transforms. 
// Reference:
//		https://www.alglib.net/fasttransforms/
library("loxxfft")
//+-------------------------------------------------------------------------------------------------------------------+
//| Fast Fourier Transform Library                                                                                      
//+-------------------------------------------------------------------------------------------------------------------+                                                                                        |
//|   void fastfouriertransform(float[] a, int nn, bool inversefft)                                                   |
//|   void realfastfouriertransform(float[] a, int tnn, bool inversefft)                                              |
//|   void fastcorrelation(float[] signal, int signallen, float[] pattern, int patternlen)                            |
//|   void fastconvolution(float[] signal, int signallen, float[] response, int negativelen, int positivelen)         |
//|   void fastsinetransform(float[] a, int tnn, bool inversefst)                                                     |
//|   void fastcosinetransform(float[] a, int tnn, bool inversefct)                                                   |
//|   void tworealffts(float[] a1, float[] a2, float[] a, float[] b, int tn)                                          |
//+-------------------------------------------------------------------------------------------------------------------+
//******************************************************************************
// Fast Fourier Transform
//
// The algorithm performs a fast Fourier transform of the complex
// function defined by nn samples on the real axis.
//
// Depending on the parameters passed, it can execute
// both direct and inverse conversion.
//
// Input parameters:
//    nn - Number of function values. Must be degree
//            deuces. Algorithm does not validate
//            passed value.
//    a - array [0 .. 2*nn-1] of Real
//            Function values. I-th value correspond
//            a[2*I] elements (real part)
//            and a[2*I+1] (the imaginary part).
//    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/
//******************************************************************************

// @function Returns Fast Fourier Transform 
// @param a float[], An array of real and imaginary parts of the function values. The real part is stored at even indices, and the imaginary part is stored at odd indices.
// @param nn int, The number of function values. It must be a power of two, but the algorithm does not validate this.
// @param inversefft bool, A boolean value that indicates the direction of the transformation. If True, it performs the inverse FFT; if False, it performs the direct FFT.
// @returns float[], Modifies the input array a in-place, which means that the transformed data (the FFT result for direct transformation or the inverse FFT result for inverse transformation) will be stored in the same array a after the function execution. The transformed data will have real and imaginary parts interleaved, with the real parts at even indices and the imaginary parts at odd indices.
export fastfouriertransform(float[] a, int nn, bool inversefft)=>
    int n = 0
    int mmax = 0
    int m = 0
    int j = 0
    int istep = 0
    int i = 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 (inversefft) 
        isign := -1
    else 
        isign := 1
 
    n := 2 * nn
    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 -= m
            m /= 2
        
        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, j - 1) - tempr)
                array.set(a, j, array.get(a, j) - tempi)
                array.set(a, i - 1, array.get(a, i) + 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)

//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
//******************************************************************************
// 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 parameters passed, it can execute
// 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/
//******************************************************************************

// @function Returns Real Fast Fourier Transform
// @param a float[], A float array containing the real-valued function samples.
// @param tnn int, The number of function values (must be a power of 2, but the algorithm does not validate this condition).
// @param inversefft bool, A boolean flag that indicates the direction of the transformation (True for inverse, False for direct).
// @returns float[], Modifies the input array a in-place, meaning that the transformed data (the FFT result for direct transformation or the inverse FFT result for inverse transformation) will be stored in the same array a after the function execution.
export 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))

//******************************************************************************
// Fast Discrete Sine Conversion
//
// The algorithm performs a fast sine transform of the real
// function defined by nn samples on the real axis.
//
// Depending on the parameters passed, it can execute
// both direct and inverse conversion.
//
// Input parameters:
//    nn - Number of function values. Must be degree
//            deuces. Algorithm does not validate
//            passed value.
//    a - array [0 .. nn-1] of Real
//            Function values.
//    Inverse FST
//        - 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/
//******************************************************************************

// @function Returns Fast Discrete Sine Conversion
// @param a float[], An array of real numbers representing the function values.
// @param tnn int, Number of function values (must be a power of two, but the code doesn't validate this).
// @param inversefst bool, A boolean flag indicating the direction of the transformation. If True, it performs the inverse FST, and if False, it performs the direct FST.
// @returns float[], The output is the transformed array 'a', which will contain the result of the transformation.
export 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)

//******************************************************************************
// Fast Discrete Cosine Transform
//
// The algorithm performs a fast cosine transform of the real
// function defined by nn samples on the real axis.
//
// Depending on the parameters passed, it can execute
// both direct and inverse conversion.
//
// Input parameters:
//    tnn - Number of function values minus one. Should be 1024
//            degree of two. The algorithm does not check
//            correct value passed.
//    a - array [0 .. nn] of Real 1025
//            Function values.
//    InverseFCT
//        - 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/
//******************************************************************************

// @function Returns Fast Discrete Cosine Transform
// @param a float[], This is a floating-point array representing the sequence of values (time-domain) that you want to transform. The function will perform the Fast Cosine Transform (FCT) or the inverse FCT on this input array, depending on the value of the inversefct parameter. The transformed result will also be stored in this same array, which means the function modifies the input array in-place.
// @param tnn int, This is an integer value representing the number of data points in the input array a. It is used to determine the size of the input array and control the loops in the algorithm. Note that the size of the input array should be a power of 2 for the Fast Cosine Transform algorithm to work correctly.
// @param inversefct bool, This is a boolean value that controls whether the function performs the regular Fast Cosine Transform or the inverse FCT. If inversefct is set to true, the function will perform the inverse FCT, and if set to false, the regular FCT will be performed. The inverse FCT can be used to transform data back into its original form (time-domain) after the regular FCT has been applied.
// @returns float[], The resulting transformed array is stored in the input array a. This means that the function modifies the input array in-place and does not return a new array.
export 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.

    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)

//******************************************************************************
// Convolution using FFT
//
// One of the folded functions is treated as a signal,
// with which we carry out convolution. The second is considered a response.
//
// At the entrance:
//    Signal - the signal with which we are convolving. array
//                    real numbers, numbering of elements
//                    0 to SignalLen-1.
//    SignalLen - signal length.
//    Response - response function. Consists of two parts,
//                    corresponding positive and negative
//                    argument values.
//
//                    Array elements with numbers from 0 to
//                    NegativeLen match response values
//                    at points from -NegativeLen to 0, respectively.
//
//                    Array elements with numbers from
//                    NegativeLen+1 to NegativeLen+PositiveLen
//                    correspond to the response values in points
//                    from 1 to PositiveLen respectively.
//                    
//    NegativeLen - "Negative length" of the response.
//    PositiveLen - "Positive length" of the response.
//                    Outside [-NegativeLen, PositiveLen]
//                    response is zero.
//                    
// At the exit:
//    Signal - function convolution values at points from 0 to
//                    SignalLen-1. https://www.alglib.net/fasttransforms/
//******************************************************************************

// @function Convolution using FFT
// @param signal float[], This is an array of real numbers representing the input signal that will be convolved with the response function. The elements are numbered from 0 to SignalLen-1.
// @param signallen int, This is an integer representing the length of the input signal array. It specifies the number of elements in the signal array.
// @param response float[], This is an array of real numbers representing the response function used for convolution. The response function consists of two parts: one corresponding to positive argument values and the other to negative argument values. Array elements with numbers from 0 to NegativeLen match the response values at points from -NegativeLen to 0, respectively. Array elements with numbers from NegativeLen+1 to NegativeLen+PositiveLen correspond to the response values in points from 1 to PositiveLen, respectively.
// @param negativelen int, This is an integer representing the "negative length" of the response function. It indicates the number of elements in the response function array that correspond to negative argument values. Outside the range [-NegativeLen, PositiveLen], the response function is considered zero.
// @param positivelen int, This is an integer representing the "positive length" of the response function. It indicates the number of elements in the response function array that correspond to positive argument values. Similar to negativelen, outside the range [-NegativeLen, PositiveLen], the response function is considered zero.
// @returns float[], The resulting convolved values are stored back in the input signal array.
export fastconvolution(float[] signal, int signallen, float[] response, int negativelen, int positivelen)=>
    int nl = 0
    int i = 0
    float t1 = 0.
    float t2 = 0.

    nl := signallen
    if negativelen > positivelen
        nl += negativelen
    else 
        nl += positivelen
    if negativelen + 1 + positivelen > nl
        nl := negativelen + 1 + positivelen
    i := 1
    while (i < nl) 
        i := i * 2
    nl := i

    float[] a1 = array.new<float>(nl, 0.)
    float[] a2 = array.new<float>(nl, 0.)

    for ix = 0 to signallen - 1 
        array.set(a1, ix, array.get(signal, ix))
    for ix = signallen to nl - 1  
        array.set(a1, ix, 0)
    for ix = 0 to nl - 1 
        array.set(a2, ix, 0)
    for ix = 1 to negativelen
        array.set(a2, nl - ix, array.get(response, negativelen - ix))

    realfastfouriertransform(a1, nl, false)
    realfastfouriertransform(a2, nl, false)

    array.set(a1, 0, array.get(a1, 0) * array.get(a2, 0))
    array.set(a1, 1, array.get(a1, 1) * array.get(a2, 1))

    for ix = 1 to nl / 2 - 1  
        t1 := array.get(a1, 2 * ix)
        t2 := array.get(a1, 2 * ix + 1)
        array.set(a1, 2 * ix, t1 * array.get(a2, 2 * ix) - t2 * array.get(a2, 2 * ix + 1))
        array.set(a1, 2 * ix + 1, t2 * array.get(a2, 2 * ix) + t1 * array.get(a2, 2 * ix + 1))
        
    realfastfouriertransform(a1, nl, true)
    
    for ix = 0 to signallen - 1 
        array.set(signal, ix, array.get(a1, ix))

//******************************************************************************
// Correlation using FFT
//
// At the entrance:
//     Signal - an array of the signal to correlate with.
//                     Item numbering from 0 to SignalLen-1
//     SignalLen - signal length.
//    
//     Pattern - an array of the pattern, the correlation of the signal with which we are looking for
//                     Numbering elements from 0 to PatternLen-1
//     PatternLen - pattern length
//
// At the exit:
//     Signal - correlation values at points from 0 to
//                     SignalLen-1. https://www.alglib.net/fasttransforms/
//******************************************************************************

// @function Returns Correlation using FFT
// @param signal float[],This is an array of real numbers representing the signal to be correlated with the pattern. The elements are numbered from 0 to SignalLen-1.
// @param signallen int, This is an integer representing the length of the input signal array.
// @param pattern float[], This is an array of real numbers representing the pattern to be correlated with the signal. The elements are numbered from 0 to PatternLen-1.
// @param patternlen int, This is an integer representing the length of the pattern array.
// @returns float[], The signal array containing the correlation values at points from 0 to SignalLen-1.
export fastcorrelation(float[] signal, int signallen, float[] pattern, int patternlen)=>
    int nl = 0
    int i = 0
    float t1 = 0.
    float t2 = 0.

    nl := signallen + patternlen
    i := 1
    while (i < nl) 
        i := i * 2
    nl := i
 
    float[] a1 = array.new<float>(nl, 0.)
    float[] a2 = array.new<float>(nl, 0.)
    
    for ix = 0 to signallen - 1 
        array.set(a1, ix, array.get(signal, ix))
    for ix = signallen to nl - 1 
        array.set(a1, ix, 0)
    for ix = 0 to patternlen - 1 
        array.set(a2, ix, array.get(pattern, ix))
    for ix = patternlen to nl - 1 
        array.set(a2, ix, 0)

    realfastfouriertransform(a1, nl, false)
    realfastfouriertransform(a2, nl, false)

    array.set(a1, 0, array.get(a1, 0) * array.get(a2, 0))
    array.set(a1, 1, array.get(a1, 1) * array.get(a2, 1))

    for ix = 1 to (nl / 2) - 1 
        t1 := array.get(a1, 2 * ix)
        t2 := array.get(a1, 2 * ix + 1)
        array.set(a1, 2 * ix, t1 * array.get(a2, 2 * ix) + t2 * array.get(a2, 2 * ix + 1))
        array.set(a1, 2 * ix + 1, t2 * array.get(a2, 2 * ix) - t1 * array.get(a2, 2 * ix + 1))

    realfastfouriertransform(a1, nl, true)
    
    for ix = 0 to signallen - 1 
        array.set(signal, ix, array.get(a1, ix))

//******************************************************************************
// Fast Fourier Transform of Two Real Functions
//
// The algorithm performs a fast Fourier transform of two
// real functions, each of which is given by tn
// readings on the real axis.
//
// The algorithm saves time, but spends only
// direct conversion.
//
// Input parameters:
//    tn - Number of function values. Must be degree
//            deuces. Algorithm does not validate
//            passed value.
//    a1 - array [0 .. nn-1] of Real
//            The values of the first function.
//    a2 - array [0 .. nn-1] of Real
//            Values of the second function.
//
// Output parameters:
//    a - Fourier transform of the first function
//    b - Fourier transform of the second function
// (for details, see the website) https://www.alglib.net/fasttransforms/
//******************************************************************************

// @function Returns Fast Fourier Transform of Two Real Functions
// @param a1 float[], An array of real numbers, representing the values of the first function.
// @param a2 float[], An array of real numbers, representing the values of the second function.
// @param a float[], An output array to store the Fourier transform of the first function.
// @param b float[], An output array to store the Fourier transform of the second function.
// @param tn float[], An integer representing the number of function values. It must be a power of two, but the algorithm doesn't validate this condition.
// @returns float[], The a and b arrays will contain the Fourier transform of the first and second functions, respectively. Note that the function overwrites the input arrays a and b.
export tworealffts(float[] a1, float[] a2, float[] a, float[] b, int tn)=>
    int j = 0
    float rep = 0.
    float rem = 0.
    float aip = 0.
    float aim = 0.
    int n = 0
    int nn = 0
    int mmax = 0
    int m = 0
    int istep = 0
    int i = 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.

    nn := tn

    for jx = 1 to tn 
        jj = jx + jx
        array.set(a, jj - 2, array.get(a1, jx - 1))
        array.set(a, jj - 1, array.get(a2, jx - 1))
    
    isign := 1
    n := 2 * nn
    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
    
    array.set(b, 0, array.get(a, 1))
    array.set(a, 1, 0.)
    array.set(b, 1, 0.)
    
    for jj = 1 to tn / 2  
        j := 2 * jj + 1
        rep := 0.5 * (array.get(a, j - 1) + array.get(a, 2 * tn + 1 - j))
        rem := 0.5 * (array.get(a, j - 1) - array.get(a, 2 * tn + 1 - j))
        aip := 0.5 * (array.get(a, j) + array.get(a, 2 * tn + 2 - j))
        aim := 0.5 * (array.get(a, j) - array.get(a, 2 * tn + 2 - j))
        array.set(a, j - 1, rep)
        array.set(a, j, aim)
        array.set(a, 2 * tn + 1 - j, rep)
        array.set(a, 2 * tn + 2 - j, -aim)
        array.set(b, j - 1, aip)
        array.set(b, j, -rem)
        array.set(b, 2 * tn + 1 - j, aip)
        array.set(b, 2 * tn + 2 - j, rem)


// example usage: RSI of Fast Discrete Cosine Transform 

float src = close
int windowper = 256 // must be power of 2
int smthcutoff = 12
int arraylevl = 0
int rsiper = 15

//check that indow is power of 2
int n = int(math.log(windowper) / math.log(2))
int N = int(math.max(math.pow(2, n), 16))

var aa = array.new<float>(N + 1, 0.)

//fill caculation array with source values
for i = 0 to N - 1
    array.set(aa, i, nz(src[i])) 

int M = array.size(aa)
int end = M - 1
smthcutoff := int(math.min(smthcutoff, M))

//regular pass of FCT
fastcosinetransform(aa, N, false)

//filter of FCT, values above smthcutoff are zeroed out
for k = 0 to end 
    if k >= smthcutoff
        array.set(aa, k, 0.)

//inverse pass of FCT
fastcosinetransform(aa, N, true)

//we are only interested in the first value; although accessing additional values can create ribbons and boundaries; end-pointed Cosine Transform
int arraylevlout = math.min(array.size(aa) - 1, arraylevl)

float out = array.get(aa, arraylevlout) 

out := ta.rsi(out, rsiper)
float sig = nz(out[1])

int mid = 50

color colorout = (out < mid)  ? color.red : (out > mid)  ? color.green : color.gray

plot(out, "RSI FDCT", color = colorout, linewidth = 2)
plot(mid, "Middle", color = bar_index % 2 ? color.gray : na)
