#_lib_FFT

//+------------------------------------------------------------------+
//|                                                       dt_FFT.mq4 |
//|                                          Copyright © 2006, klot. |
//|                                                     klot@mail.ru |
//+------------------------------------------------------------------+
#property copyright "FFT for MQL4 klot"
#property link      "http://alglib.sources.ru/fft/"
#property library

#define pi 3.14159265358979323846
//+-------------------------------------------------------------------------------------------------------------------+
//| My function                                                                                                       |
//+-------------------------------------------------------------------------------------------------------------------+
//|   #import "#_lib_FFT.ex4"                                                                                         |
//|   void fastfouriertransform(double& a[], int nn, bool inversefft);                                                |
//|   void realfastfouriertransform(double& a[], int tnn, bool inversefft);                                           |
//|   void fastcorellation(double& signal[], int signallen,  double& pattern[],  int patternlen);                     |
//|   void fastconvolution(double& signal[], int signallen, double& response[], int negativelen, int positivelen);    |
//|   void fastsinetransform(double& a[], int tnn, bool inversefst);                                                  |
//|   void fastcosinetransform(double& a[], int tnn, bool inversefct);                                                |
//|   void tworealffts(double a1[], double a2[],  double& a[],  double& b[], int tn);                                 |
//|   #import                                                                                                         |
//+-------------------------------------------------------------------------------------------------------------------+
/*************************************************************************
Áûñòðîå ïðåîáðàçîâàíèå Ôóðüå

Àëãîðèòì ïðîâîäèò áûñòðîå ïðåîáðàçîâàíèå Ôóðüå êîìïëåêñíîé
ôóíêöèè, çàäàííîé nn îòñ÷åòàìè íà äåéñòâèòåëüíîé îñè.

 çàâèñèìîñòè îò  ïåðåäàííûõ ïàðàìåòðîâ, ìîæåò âûïîëíÿòüñÿ
êàê ïðÿìîå, òàê è îáðàòíîå ïðåîáðàçîâàíèå.

Âõîäíûå ïàðàìåòðû:
    nn  -   ×èñëî çíà÷åíèé ôóíêöèè. Äîëæíî  áûòü  ñòåïåíüþ
            äâîéêè. Àëãîðèòì   íå  ïðîâåðÿåò  ïðàâèëüíîñòü
            ïåðåäàííîãî çíà÷åíèÿ.
    a   -   array [0 .. 2*nn-1] of Real
            Çíà÷åíèÿ ôóíêöèè. I-îìó çíà÷åíèþ ñîîòâåòñòâóþò
            ýëåìåíòû   a[2*I]     (âåùåñòâåííàÿ     ÷àñòü)
            è a[2*I+1] (ìíèìàÿ ÷àñòü).
    InverseFFT
        -   íàïðàâëåíèå ïðåîáðàçîâàíèÿ.
            True, åñëè îáðàòíîå, False, åñëè ïðÿìîå.
            
Âûõîäíûå ïàðàìåòðû:
    a   -   ðåçóëüòàò   ïðåîáðàçîâàíèÿ.   Ïîäðîáíåå    ñì.
            îïèñàíèå íà ñàéòå.  http://alglib.sources.ru/fft/
*************************************************************************/
void fastfouriertransform(double& a[], int nn, bool inversefft)
{
    int ii;
    int jj;
    int n;
    int mmax;
    int m;
    int j;
    int istep;
    int i;
    int isign;
    double wtemp;
    double wr;
    double wpr;
    double wpi;
    double wi;
    double theta;
    double tempr;
    double tempi;

    if( inversefft )
    {
        isign = -1;
    }
    else
    {
        isign = 1;
    }
    n = 2*nn;
    j = 1;
    for(ii = 1; ii <= nn; ii++)
    {
        i = 2*ii-1;
        if( j>i )
        {
            tempr = a[j-1];
            tempi = a[j];
            a[j-1] = a[i-1];
            a[j] = a[i];
            a[i-1] = tempr;
            a[i] = tempi;
        }
        m = n/2;
        while(m>=2&&j>m)
        {
            j = j-m;
            m = m/2;
        }
        j = j+m;
    }
    mmax = 2;
    while(n>mmax)
    {
        istep = 2*mmax;
        theta = 2.0*pi/(isign*mmax);
        wpr = -2.0*MathPow(MathSin(0.5*theta),2);
        wpi = MathSin(theta);
        wr = 1.0;
        wi = 0.0;
        for(ii = 1; ii <= mmax/2; ii++)
        {
            m = 2*ii-1;
            for(jj = 0; jj <= (n-m)/istep; jj++)
            {
                i = m+jj*istep;
                j = i+mmax;
                tempr = wr*a[j-1]-wi*a[j];
                tempi = wr*a[j]+wi*a[j-1];
                a[j-1] = a[i-1]-tempr;
                a[j] = a[i]-tempi;
                a[i-1] = a[i-1]+tempr;
                a[i] = a[i]+tempi;
            }
            wtemp = wr;
            wr = wr*wpr-wi*wpi+wr;
            wi = wi*wpr+wtemp*wpi+wi;
        }
        mmax = istep;
    }
    if( inversefft )
    {
        for(i = 1; i <= 2*nn; i++)
        {
            a[i-1] = a[i-1]/nn;
        }
    }
    return;
}
/*************************************************************************
Áûñòðîå ïðåîáðàçîâàíèå Ôóðüå

Àëãîðèòì ïðîâîäèò áûñòðîå ïðåîáðàçîâàíèå Ôóðüå âåùåñòâåííîé
ôóíêöèè, çàäàííîé n îòñ÷åòàìè íà äåéñòâèòåëüíîé îñè.

 çàâèñèìîñòè îò  ïåðåäàííûõ ïàðàìåòðîâ, ìîæåò âûïîëíÿòüñÿ
êàê ïðÿìîå, òàê è îáðàòíîå ïðåîáðàçîâàíèå.

Âõîäíûå ïàðàìåòðû:
    tnn  -   ×èñëî çíà÷åíèé ôóíêöèè. Äîëæíî  áûòü  ñòåïåíüþ
            äâîéêè. Àëãîðèòì   íå  ïðîâåðÿåò  ïðàâèëüíîñòü
            ïåðåäàííîãî çíà÷åíèÿ.
    a   -   array [0 .. nn-1] of Real
            Çíà÷åíèÿ ôóíêöèè.
    InverseFFT
        -   íàïðàâëåíèå ïðåîáðàçîâàíèÿ.
            True, åñëè îáðàòíîå, False, åñëè ïðÿìîå.
            
Âûõîäíûå ïàðàìåòðû:
    a   -   ðåçóëüòàò   ïðåîáðàçîâàíèÿ.   Ïîäðîáíåå    ñì.
            îïèñàíèå íà ñàéòå. http://alglib.sources.ru/fft/
*************************************************************************/
void realfastfouriertransform(double& a[], int tnn, bool inversefft)
{
    double twr;
    double twi;
    double twpr;
    double twpi;
    double twtemp;
    double ttheta;
    int i;
    int i1;
    int i2;
    int i3;
    int i4;
    double c1;
    double c2;
    double h1r;
    double h1i;
    double h2r;
    double h2i;
    double wrs;
    double wis;
    int nn;
    int ii;
    int jj;
    int n;
    int mmax;
    int m;
    int j;
    int istep;
    int isign;
    double wtemp;
    double wr;
    double wpr;
    double wpi;
    double wi;
    double theta;
    double tempr;
    double tempi;

    if( tnn==1 )
    {
        return;
    }
    if( !inversefft )
    {
        ttheta = 2.0*pi/tnn;
        c1 = 0.5;
        c2 = -0.5;
    }
    else
    {
        ttheta = 2.0*pi/tnn;
        c1 = 0.5;
        c2 = 0.5;
        ttheta = -ttheta;
        twpr = -2.0*MathPow(MathSin(0.5*ttheta),2);
        twpi = MathSin(ttheta);
        twr = 1.0+twpr;
        twi = twpi;
        for(i = 2; i <= tnn/4+1; i++)
        {
            i1 = i+i-2;
            i2 = i1+1;
            i3 = tnn+1-i2;
            i4 = i3+1;
            wrs = twr;
            wis = twi;
            h1r = c1*(a[i1]+a[i3]);
            h1i = c1*(a[i2]-a[i4]);
            h2r = -c2*(a[i2]+a[i4]);
            h2i = c2*(a[i1]-a[i3]);
            a[i1] = h1r+wrs*h2r-wis*h2i;
            a[i2] = h1i+wrs*h2i+wis*h2r;
            a[i3] = h1r-wrs*h2r+wis*h2i;
            a[i4] = -h1i+wrs*h2i+wis*h2r;
            twtemp = twr;
            twr = twr*twpr-twi*twpi+twr;
            twi = twi*twpr+twtemp*twpi+twi;
        }
        h1r = a[0];
        a[0] = c1*(h1r+a[1]);
        a[1] = c1*(h1r-a[1]);
    }
    if( inversefft )
    {
        isign = -1;
    }
    else
    {
        isign = 1;
    }
    n = tnn;
    nn = tnn/2;
    j = 1;
    for(ii = 1; ii <= nn; ii++)
    {
        i = 2*ii-1;
        if( j>i )
        {
            tempr = a[j-1];
            tempi = a[j];
            a[j-1] = a[i-1];
            a[j] = a[i];
            a[i-1] = tempr;
            a[i] = tempi;
        }
        m = n/2;
        while(m>=2&&j>m)
        {
            j = j-m;
            m = m/2;
        }
        j = j+m;
    }
    mmax = 2;
    while(n>mmax)
    {
        istep = 2*mmax;
        theta = 2.0*pi/(isign*mmax);
        wpr = -2.0*MathPow(MathSin(0.5*theta),2);
        wpi = MathSin(theta);
        wr = 1.0;
        wi = 0.0;
        for(ii = 1; ii <= mmax/2; ii++)
        {
            m = 2*ii-1;
            for(jj = 0; jj <= (n-m)/istep; jj++)
            {
                i = m+jj*istep;
                j = i+mmax;
                tempr = wr*a[j-1]-wi*a[j];
                tempi = wr*a[j]+wi*a[j-1];
                a[j-1] = a[i-1]-tempr;
                a[j] = a[i]-tempi;
                a[i-1] = a[i-1]+tempr;
                a[i] = a[i]+tempi;
            }
            wtemp = wr;
            wr = wr*wpr-wi*wpi+wr;
            wi = wi*wpr+wtemp*wpi+wi;
        }
        mmax = istep;
    }
    if( inversefft )
    {
        for(i = 1; i <= 2*nn; i++)
        {
            a[i-1] = a[i-1]/nn;
        }
    }
    if( !inversefft )
    {
        twpr = -2.0*MathPow(MathSin(0.5*ttheta),2);
        twpi = MathSin(ttheta);
        twr = 1.0+twpr;
        twi = twpi;
        for(i = 2; i <= tnn/4+1; i++)
        {
            i1 = i+i-2;
            i2 = i1+1;
            i3 = tnn+1-i2;
            i4 = i3+1;
            wrs = twr;
            wis = twi;
            h1r = c1*(a[i1]+a[i3]);
            h1i = c1*(a[i2]-a[i4]);
            h2r = -c2*(a[i2]+a[i4]);
            h2i = c2*(a[i1]-a[i3]);
            a[i1] = h1r+wrs*h2r-wis*h2i;
            a[i2] = h1i+wrs*h2i+wis*h2r;
            a[i3] = h1r-wrs*h2r+wis*h2i;
            a[i4] = -h1i+wrs*h2i+wis*h2r;
            twtemp = twr;
            twr = twr*twpr-twi*twpi+twr;
            twi = twi*twpr+twtemp*twpi+twi;
        }
        h1r = a[0];
        a[0] = h1r+a[1];
        a[1] = h1r-a[1];
    }
    return;
}
//===========================================================================================
/*************************************************************************
Êîðåëëÿöèÿ ñ èñïîëüçîâàíèåì ÁÏÔ

Íà âõîäå:
    Signal      -   ìàññèâ ñèãíàë, ñ êîòîðûì ïðîâîäèì êîðåëëÿöèþ.
                    Íóìåðàöèÿ ýëåìåíòîâ îò 0 äî SignalLen-1
    SignalLen   -   äëèíà ñèãíàëà.
    
    Pattern     -   ìàññèâ îáðàçåö, êîðåëëÿöèþ ñèãíàëà ñ êîòîðûì ìû èùåì
                    Íóìåðàöèÿ ýëåìåíòîâ îò 0 äî PatternLen-1
    PatternLen  -   äëèíà îáðàçöà

Íà âûõîäå:
    Signal      -   çíà÷åíèÿ êîðåëëÿöèè â òî÷êàõ îò 0 äî
                    SignalLen-1. http://alglib.sources.ru/fft/
*************************************************************************/
void fastcorellation(double& signal[], int signallen,  double pattern[],  int patternlen)
{
    //ap::real_1d_array a1;
    //ap::real_1d_array a2;
    double a1[], a2[];
    //---
    int nl;
    int i;
    double t1;
    double t2;

    nl = signallen+patternlen;
    i = 1;
    while(i<nl)
    {
        i = i*2;
    }
    nl = i;
    //a1.setbounds(0, nl-1);
    //a2.setbounds(0, nl-1);
    ArrayResize(a1, nl);
    ArrayResize(a2, nl);
    //---
    for(i = 0; i <= signallen-1; i++)
    {
        a1[i] = signal[i];
    }
    for(i = signallen; i <= nl-1; i++)
    {
        a1[i] = 0;
    }
    for(i = 0; i <= patternlen-1; i++)
    {
        a2[i] = pattern[i];
    }
    for(i = patternlen; i <= nl-1; i++)
    {
        a2[i] = 0;
    }
    realfastfouriertransform(a1, nl, false);
    realfastfouriertransform(a2, nl, false);
    a1[0] = a1[0]*a2[0];
    a1[1] = a1[1]*a2[1];
    for(i = 1; i <= (nl/2)-1; i++)
    {
        t1 = a1[2*i];
        t2 = a1[2*i+1];
        a1[2*i] = t1*a2[2*i]+t2*a2[2*i+1];
        a1[2*i+1] = t2*a2[2*i]-t1*a2[2*i+1];
    }
    realfastfouriertransform(a1, nl, true);
    for(i = 0; i <= signallen-1; i++)
    {
        signal[i] = a1[i];
    }
    return;
}
//==========================================================================================
/*************************************************************************
Ñâåðòêà ñ èñïîëüçîâàíèåì ÁÏÔ

Îäíà èç ñâîðà÷èâàåìûõ ôóíêöèé òðàêòóåòñÿ, êàê ñèãíàë,
ñ êîòîðûì ïðîâîäèì ñâåðòêó. Âòîðàÿ ñ÷èòàåòñÿ îòêëèêîì.

Íà âõîäå:
    Signal      -   ñèãíàë, ñ êîòîðûì ïðîâîäèì ñâåðòêó. Ìàññèâ
                    âåùåñòâåííûõ  ÷èñåë,  íóìåðàöèÿ  ýëåìåíòîâ
                    îò 0 äî SignalLen-1.
    SignalLen   -   äëèíà ñèãíàëà.
    Response    -   ôóíêöèÿ îòêëèêà. Ñîñòîèò èç  äâóõ  ÷àñòåé,
                    ñîîòâåòñòâóþùèõ ïîëîæèòåëüíûì è îòðèöàòåëüíûì
                    çíà÷åíèÿì àðãóìåíòà.

                    Ýëåìåíòàì    ìàññèâà  ñ   íîìåðàìè îò 0 äî
                    NegativeLen ñîîòâåòñòâóþò çíà÷åíèÿ îòêëèêà
                    â òî÷êàõ îò -NegativeLen äî 0 ñîîòâåòñòâåííî.

                    Ýëåìåíòàì      ìàññèâà   ñ   íîìåðàìè   îò
                    NegativeLen+1  äî  NegativeLen+PositiveLen
                    ñîîòâåòñòâóþò  çíà÷åíèÿ  îòêëèêà â  òî÷êàõ
                    îò 1 äî PositiveLen ñîîòâåòñòâåííî.
                    
    NegativeLen -   "Îòðèöàòåëüíàÿ äëèíà" îòêëèêà.
    PositiveLen -   "Ïîëîæèòåëüíàÿ äëèíà" îòêëèêà.
                    Çà   ïðåäåëàìè [-NegativeLen, PositiveLen]
                    îòêëèê ðàâåí íîëþ.
                    
Íà âûõîäå:
    Signal      -   çíà÷åíèÿ ñâåðòêè ôóíêöèè â òî÷êàõ îò 0  äî
                    SignalLen-1. http://alglib.sources.ru/fft/
*************************************************************************/
void fastconvolution(double& signal[], int signallen, double& response[], int negativelen, int positivelen)
{
    //ap::real_1d_array a1;
    //ap::real_1d_array a2;
    double a1[], a2[];
    //----
    int nl;
    int i;
    double t1;
    double t2;

    nl = signallen;
    if( negativelen>positivelen )
    {
        nl = nl+negativelen;
    }
    else
    {
        nl = nl+positivelen;
    }
    if( negativelen+1+positivelen>nl )
    {
        nl = negativelen+1+positivelen;
    }
    i = 1;
    while(i<nl)
    {
        i = i*2;
    }
    nl = i;
    //a1.setbounds(0, nl-1);
    //a2.setbounds(0, nl-1);
    ArrayResize(a1, nl);
    ArrayResize(a2, nl);
    //---
    for(i = 0; i <= signallen-1; i++)
    {
        a1[i] = signal[i];
    }
    for(i = signallen; i <= nl-1; i++)
    {
        a1[i] = 0;
    }
    for(i = 0; i <= nl-1; i++)
    {
        a2[i] = 0;
    }
    for(i = 0; i <= positivelen; i++)
    {
        a2[i] = response[i+negativelen];
    }
    for(i = 1; i <= negativelen; i++)
    {
        a2[nl-i] = response[negativelen-i];
    }
    realfastfouriertransform(a1, nl, false);
    realfastfouriertransform(a2, nl, false);
    a1[0] = a1[0]*a2[0];
    a1[1] = a1[1]*a2[1];
    for(i = 1; i <= nl/2-1; i++)
    {
        t1 = a1[2*i];
        t2 = a1[2*i+1];
        a1[2*i] = t1*a2[2*i]-t2*a2[2*i+1];
        a1[2*i+1] = t2*a2[2*i]+t1*a2[2*i+1];
    }
    realfastfouriertransform(a1, nl, true);
    for(i = 0; i <= signallen-1; i++)
    {
        signal[i] = a1[i];
    }
    return;
}
//===========================================================================================
/*************************************************************************
Áûñòðîå äèñêðåòíîå ñèíóñíîå ïðåîáðàçîâàíèå

Àëãîðèòì ïðîâîäèò áûñòðîå ñèíóñíîå ïðåîáðàçîâàíèå âåùåñòâåííîé
ôóíêöèè, çàäàííîé nn îòñ÷åòàìè íà äåéñòâèòåëüíîé îñè.

 çàâèñèìîñòè îò  ïåðåäàííûõ ïàðàìåòðîâ, ìîæåò âûïîëíÿòüñÿ
êàê ïðÿìîå, òàê è îáðàòíîå ïðåîáðàçîâàíèå.

Âõîäíûå ïàðàìåòðû:
    nn  -   ×èñëî çíà÷åíèé ôóíêöèè. Äîëæíî  áûòü  ñòåïåíüþ
            äâîéêè. Àëãîðèòì   íå  ïðîâåðÿåò  ïðàâèëüíîñòü
            ïåðåäàííîãî çíà÷åíèÿ.
    a   -   array [0 .. nn-1] of Real
            Çíà÷åíèÿ ôóíêöèè.
    InverseFST
        -   íàïðàâëåíèå ïðåîáðàçîâàíèÿ.
            True, åñëè îáðàòíîå, False, åñëè ïðÿìîå.
            
Âûõîäíûå ïàðàìåòðû:
    a   -   ðåçóëüòàò   ïðåîáðàçîâàíèÿ.   Ïîäðîáíåå    ñì.
            îïèñàíèå íà ñàéòå. http://alglib.sources.ru/fft/
*************************************************************************/
void fastsinetransform(double& a[], int tnn, bool inversefst)
{
    int jj;
    int j;
    int tm;
    int n2;
    double sum;
    double y1;
    double y2;
    double theta;
    double wi;
    double wr;
    double wpi;
    double wpr;
    double wtemp;
    double twr;
    double twi;
    double twpr;
    double twpi;
    double twtemp;
    double ttheta;
    int i;
    int i1;
    int i2;
    int i3;
    int i4;
    double c1;
    double c2;
    double h1r;
    double h1i;
    double h2r;
    double h2i;
    double wrs;
    double wis;
    int nn;
    int ii;
    int n;
    int mmax;
    int m;
    int istep;
    int isign;
    double tempr;
    double tempi;

    if( tnn==1 )
    {
        a[0] = 0;
        return;
    }
    theta = pi/tnn;
    wr = 1.0;
    wi = 0.0;
    wpr = -2.0*MathPow(MathSin(0.5*theta),2);
    wpi = MathSin(theta);
    a[0] = 0.0;
    tm = tnn/2;
    n2 = tnn+2;
    for(j = 2; j <= tm+1; j++)
    {
        wtemp = wr;
        wr = wr*wpr-wi*wpi+wr;
        wi = wi*wpr+wtemp*wpi+wi;
        y1 = wi*(a[j-1]+a[n2-j-1]);
        y2 = 0.5*(a[j-1]-a[n2-j-1]);
        a[j-1] = y1+y2;
        a[n2-j-1] = y1-y2;
    }
    ttheta = 2.0*pi/tnn;
    c1 = 0.5;
    c2 = -0.5;
    isign = 1;
    n = tnn;
    nn = tnn/2;
    j = 1;
    for(ii = 1; ii <= nn; ii++)
    {
        i = 2*ii-1;
        if( j>i )
        {
            tempr = a[j-1];
            tempi = a[j];
            a[j-1] = a[i-1];
            a[j] = a[i];
            a[i-1] = tempr;
            a[i] = tempi;
        }
        m = n/2;
        while(m>=2&&j>m)
        {
            j = j-m;
            m = m/2;
        }
        j = j+m;
    }
    mmax = 2;
    while(n>mmax)
    {
        istep = 2*mmax;
        theta = 2.0*pi/(isign*mmax);
        wpr = -2.0*MathPow(MathSin(0.5*theta),2);
        wpi = MathSin(theta);
        wr = 1.0;
        wi = 0.0;
        for(ii = 1; ii <= mmax/2; ii++)
        {
            m = 2*ii-1;
            for(jj = 0; jj <= (n-m)/istep; jj++)
            {
                i = m+jj*istep;
                j = i+mmax;
                tempr = wr*a[j-1]-wi*a[j];
                tempi = wr*a[j]+wi*a[j-1];
                a[j-1] = a[i-1]-tempr;
                a[j] = a[i]-tempi;
                a[i-1] = a[i-1]+tempr;
                a[i] = a[i]+tempi;
            }
            wtemp = wr;
            wr = wr*wpr-wi*wpi+wr;
            wi = wi*wpr+wtemp*wpi+wi;
        }
        mmax = istep;
    }
    twpr = -2.0*MathPow(MathSin(0.5*ttheta),2);
    twpi = MathSin(ttheta);
    twr = 1.0+twpr;
    twi = twpi;
    for(i = 2; i <= tnn/4+1; i++)
    {
        i1 = i+i-2;
        i2 = i1+1;
        i3 = tnn+1-i2;
        i4 = i3+1;
        wrs = twr;
        wis = twi;
        h1r = c1*(a[i1]+a[i3]);
        h1i = c1*(a[i2]-a[i4]);
        h2r = -c2*(a[i2]+a[i4]);
        h2i = c2*(a[i1]-a[i3]);
        a[i1] = h1r+wrs*h2r-wis*h2i;
        a[i2] = h1i+wrs*h2i+wis*h2r;
        a[i3] = h1r-wrs*h2r+wis*h2i;
        a[i4] = -h1i+wrs*h2i+wis*h2r;
        twtemp = twr;
        twr = twr*twpr-twi*twpi+twr;
        twi = twi*twpr+twtemp*twpi+twi;
    }
    h1r = a[0];
    a[0] = h1r+a[1];
    a[1] = h1r-a[1];
    sum = 0.0;
    a[0] = 0.5*a[0];
    a[1] = 0.0;
    for(jj = 0; jj <= tm-1; jj++)
    {
        j = 2*jj+1;
        sum = sum+a[j-1];
        a[j-1] = a[j];
        a[j] = sum;
    }
    if( inversefst )
    {
        for(j = 1; j <= tnn; j++)
        {
            a[j-1] = a[j-1]*2/tnn;
        }
    }
    return;
}
//=============================================================================================
/*************************************************************************
Áûñòðîå äèñêðåòíîå êîñèíóñíîå ïðåîáðàçîâàíèå

Àëãîðèòì ïðîâîäèò áûñòðîå êîñèíóñíîå ïðåîáðàçîâàíèå âåùåñòâåííîé
ôóíêöèè, çàäàííîé nn îòñ÷åòàìè íà äåéñòâèòåëüíîé îñè.

 çàâèñèìîñòè îò  ïåðåäàííûõ ïàðàìåòðîâ, ìîæåò âûïîëíÿòüñÿ
êàê ïðÿìîå, òàê è îáðàòíîå ïðåîáðàçîâàíèå.

Âõîäíûå ïàðàìåòðû:
    tnn  -  ×èñëî çíà÷åíèé ôóíêöèè ìèíóñ îäèí. Äîëæíî áûòü                 1024
            ñòåïåíüþ   äâîéêè.   Àëãîðèòì   íå   ïðîâåðÿåò
            ïðàâèëüíîñòü ïåðåäàííîãî çíà÷åíèÿ.
    a   -   array [0 .. nn] of Real                                        1025
            Çíà÷åíèÿ ôóíêöèè.
    InverseFCT
        -   íàïðàâëåíèå ïðåîáðàçîâàíèÿ.
            True, åñëè îáðàòíîå, False, åñëè ïðÿìîå.
            
Âûõîäíûå ïàðàìåòðû:
    a   -   ðåçóëüòàò   ïðåîáðàçîâàíèÿ.   Ïîäðîáíåå    ñì.
            îïèñàíèå íà ñàéòå. http://alglib.sources.ru/fft/
*************************************************************************/
void fastcosinetransform(double& a[], int tnn, bool inversefct)
{
    int j;
    int n2;
    double sum;
    double y1;
    double y2;
    double theta;
    double wi;
    double wpi;
    double wr;
    double wpr;
    double wtemp;
    double twr;
    double twi;
    double twpr;
    double twpi;
    double twtemp;
    double ttheta;
    int i;
    int i1;
    int i2;
    int i3;
    int i4;
    double c1;
    double c2;
    double h1r;
    double h1i;
    double h2r;
    double h2i;
    double wrs;
    double wis;
    int nn;
    int ii;
    int jj;
    int n;
    int mmax;
    int m;
    int istep;
    int isign;
    double tempr;
    double tempi;

    if( tnn==1 )
    {
        y1 = a[0];
        y2 = a[1];
        a[0] = 0.5*(y1+y2);
        a[1] = 0.5*(y1-y2);
        if( inversefct )
        {
            a[0] = a[0]*2;
            a[1] = a[1]*2;
        }
        return;
    }
    wi = 0;
    wr = 1;
    theta = pi/tnn;
    wtemp = MathSin(theta*0.5);
    wpr = -2.0*wtemp*wtemp;
    wpi = MathSin(theta);
    sum = 0.5*(a[0]-a[tnn]);
    a[0] = 0.5*(a[0]+a[tnn]);
    n2 = tnn+2;
    for(j = 2; j <= tnn/2; j++)
    {
        wtemp = wr;
        wr = wtemp*wpr-wi*wpi+wtemp;
        wi = wi*wpr+wtemp*wpi+wi;
        y1 = 0.5*(a[j-1]+a[n2-j-1]);
        y2 = a[j-1]-a[n2-j-1];
        a[j-1] = y1-wi*y2;
        a[n2-j-1] = y1+wi*y2;
        sum = sum+wr*y2;
    }
    ttheta = 2.0*pi/tnn;
    c1 = 0.5;
    c2 = -0.5;
    isign = 1;
    n = tnn;
    nn = tnn/2;
    j = 1;
    for(ii = 1; ii <= nn; ii++)
    {
        i = 2*ii-1;
        if( j>i )
        {
            tempr = a[j-1];
            tempi = a[j];
            a[j-1] = a[i-1];
            a[j] = a[i];
            a[i-1] = tempr;
            a[i] = tempi;
        }
        m = n/2;
        while(m>=2&&j>m)
        {
            j = j-m;
            m = m/2;
        }
        j = j+m;
    }
    mmax = 2;
    while(n>mmax)
    {
        istep = 2*mmax;
        theta = 2.0*pi/(isign*mmax);
        wpr = -2.0*MathPow(MathSin(0.5*theta),2);
        wpi = MathSin(theta);
        wr = 1.0;
        wi = 0.0;
        for(ii = 1; ii <= mmax/2; ii++)
        {
            m = 2*ii-1;
            for(jj = 0; jj <= (n-m)/istep; jj++)
            {
                i = m+jj*istep;
                j = i+mmax;
                tempr = wr*a[j-1]-wi*a[j];
                tempi = wr*a[j]+wi*a[j-1];
                a[j-1] = a[i-1]-tempr;
                a[j] = a[i]-tempi;
                a[i-1] = a[i-1]+tempr;
                a[i] = a[i]+tempi;
            }
            wtemp = wr;
            wr = wr*wpr-wi*wpi+wr;
            wi = wi*wpr+wtemp*wpi+wi;
        }
        mmax = istep;
    }
    twpr = -2.0*MathPow(MathSin(0.5*ttheta),2);
    twpi = MathSin(ttheta);
    twr = 1.0+twpr;
    twi = twpi;
    for(i = 2; i <= tnn/4+1; i++)
    {
        i1 = i+i-2;
        i2 = i1+1;
        i3 = tnn+1-i2;
        i4 = i3+1;
        wrs = twr;
        wis = twi;
        h1r = c1*(a[i1]+a[i3]);
        h1i = c1*(a[i2]-a[i4]);
        h2r = -c2*(a[i2]+a[i4]);
        h2i = c2*(a[i1]-a[i3]);
        a[i1] = h1r+wrs*h2r-wis*h2i;
        a[i2] = h1i+wrs*h2i+wis*h2r;
        a[i3] = h1r-wrs*h2r+wis*h2i;
        a[i4] = -h1i+wrs*h2i+wis*h2r;
        twtemp = twr;
        twr = twr*twpr-twi*twpi+twr;
        twi = twi*twpr+twtemp*twpi+twi;
    }
    h1r = a[0];
    a[0] = h1r+a[1];
    a[1] = h1r-a[1];
    a[tnn] = a[1];
    a[1] = sum;
    j = 4;
    while(j<=tnn)
    {
        sum = sum+a[j-1];
        a[j-1] = sum;
        j = j+2;
    }
    if( inversefct )
    {
        for(j = 0; j <= tnn; j++)
        {
            a[j] = a[j]*2/tnn;
        }
    }
    return;
}
//===============================================================================================
/*************************************************************************
Áûñòðîå ïðåîáðàçîâàíèå Ôóðüå äâóõ âåùåñòâåííûõ ôóíêöèé

Àëãîðèòì ïðîâîäèò   áûñòðîå   ïðåîáðàçîâàíèå   Ôóðüå  äâóõ
âåùåñòâåííûõ    ôóíêöèé,  êàæäàÿ   èç  êîòîðûõ  çàäàíà  tn
îòñ÷åòàìè íà äåéñòâèòåëüíîé îñè.

Àëãîðèòì  ïîçâîëÿåò  ñýêîíîìèòü  âðåìÿ, íî ïðîâîäèò òîëüêî
ïðÿìîå ïðåîáðàçîâàíèå.

Âõîäíûå ïàðàìåòðû:
    tn  -   ×èñëî çíà÷åíèé ôóíêöèé. Äîëæíî  áûòü  ñòåïåíüþ
            äâîéêè. Àëãîðèòì   íå  ïðîâåðÿåò  ïðàâèëüíîñòü
            ïåðåäàííîãî çíà÷åíèÿ.
    a1  -   array [0 .. nn-1] of Real
            Çíà÷åíèÿ ïåðâîé ôóíêöèè.
    a2  -   array [0 .. nn-1] of Real
            Çíà÷åíèÿ âòîðîé ôóíêöèè.

Âûõîäíûå ïàðàìåòðû:
    a   -   Ïðåîáðàçîâàíèå Ôóðüå ïåðâîé ôóíêöèè
    b   -   Ïðåîáðàçîâàíèå Ôóðüå âòîðîé ôóíêöèè
(ïîäðîáíåå ñì. íà ñàéòå) http://alglib.sources.ru/fft/
*************************************************************************/
void tworealffts(double a1[], double a2[],  double& a[],  double& b[], int tn)
{
    int jj;
    int j;
    double rep;
    double rem;
    double aip;
    double aim;
    int ii;
    int n;
    int nn;
    int mmax;
    int m;
    int istep;
    int i;
    int isign;
    double wtemp;
    double wr;
    double wpr;
    double wpi;
    double wi;
    double theta;
    double tempr;
    double tempi;

    nn = tn;
    //a.setbounds(0, 2*tn-1);
    //b.setbounds(0, 2*tn-1);
    ArrayResize(a, 2*tn);
    ArrayResize(b, 2*tn);
    //---
    for(j = 1; j <= tn; j++)
    {
        jj = j+j;
        a[jj-2] = a1[j-1];
        a[jj-1] = a2[j-1];
    }
    isign = 1;
    n = 2*nn;
    j = 1;
    for(ii = 1; ii <= nn; ii++)
    {
        i = 2*ii-1;
        if( j>i )
        {
            tempr = a[j-1];
            tempi = a[j];
            a[j-1] = a[i-1];
            a[j] = a[i];
            a[i-1] = tempr;
            a[i] = tempi;
        }
        m = n/2;
        while(m>=2&&j>m)
        {
            j = j-m;
            m = m/2;
        }
        j = j+m;
    }
    mmax = 2;
    while(n>mmax)
    {
        istep = 2*mmax;
        theta = 2.0*pi/(isign*mmax);
        wpr = -2.0*MathPow(MathSin(0.5*theta),2);
        wpi = MathSin(theta);
        wr = 1.0;
        wi = 0.0;
        for(ii = 1; ii <= mmax/2; ii++)
        {
            m = 2*ii-1;
            for(jj = 0; jj <= (n-m)/istep; jj++)
            {
                i = m+jj*istep;
                j = i+mmax;
                tempr = wr*a[j-1]-wi*a[j];
                tempi = wr*a[j]+wi*a[j-1];
                a[j-1] = a[i-1]-tempr;
                a[j] = a[i]-tempi;
                a[i-1] = a[i-1]+tempr;
                a[i] = a[i]+tempi;
            }
            wtemp = wr;
            wr = wr*wpr-wi*wpi+wr;
            wi = wi*wpr+wtemp*wpi+wi;
        }
        mmax = istep;
    }
    b[0] = a[1];
    a[1] = 0.0;
    b[1] = 0.0;
    for(jj = 1; jj <= tn/2; jj++)
    {
        j = 2*jj+1;
        rep = 0.5*(a[j-1]+a[2*tn+1-j]);
        rem = 0.5*(a[j-1]-a[2*tn+1-j]);
        aip = 0.5*(a[j]+a[2*tn+2-j]);
        aim = 0.5*(a[j]-a[2*tn+2-j]);
        a[j-1] = rep;
        a[j] = aim;
        a[2*tn+1-j] = rep;
        a[2*tn+2-j] = -aim;
        b[j-1] = aip;
        b[j] = -rem;
        b[2*tn+1-j] = aip;
        b[2*tn+2-j] = rem;
    }
    return;
}
//==============================================================================================