/*====================================================================================
    EVS Codec 3GPP TS26.443 Jun 30, 2015. Version CR 26.443-0006
  ====================================================================================*/

#include <math.h>
#include "cnst.h"
#include "options.h"
#include "prot.h"

/*-------------------------------------------------------------------*
 * four1()
 *
 *  From "numerical recipes in C".
 *  Replace data by its DFT, if isign is input as 1; or replace data
 *  by nn times its inverse-DFT, if isign is input as -1.
 *  data is a complex array of length nn, input as a real
 *  array data[1...2nn]. nn must be an integer power of 2
 *-------------------------------------------------------------------*/

static void four1(
    float *data,     /* i/o: data array   .......... */
    short nn,        /* i  : length of data array    */
    short isign      /* i  : sign +1 or -1           */
)
{
    short n,mmax,m,j,istep,i;
    float wtemp,wr,wpr,wpi,wi,theta;
    float tempr,tempi;

    n=nn << 1;
    j=1;

    /* this is the bit-reversal section of the routine */
    for (i=1; i<n; i+=2)
    {
        if (j > i)
        {
            /* exchange the two complex numbers */
            SWAP(data[j],data[i]);
            SWAP(data[j+1],data[i+1]);
        }
        m=n >> 1;
        while (m >= 2 && j > m)
        {
            j -= m;
            m >>= 1;
        }
        j += m;
    }
    mmax=2;
    /* here begins the Danielson-Lanczos section of the routine */
    /* Outer loop executed log2(nn) times */
    while (n > mmax)
    {
        istep=2*mmax;
        /* initialization for the trigonometric recurrence */
        theta=(float) (6.28318530717959/(isign*mmax));
        wtemp=(float) (sin(0.5f*theta));
        wpr = -2.0f*wtemp*wtemp;
        wpi=(float) sin(theta);
        wr=1.0f;
        wi=0.0f;
        /* here are the two nested loops */
        for (m=1; m<mmax; m+=2)
        {
            for (i=m; i<=n; i+=istep)
            {
                /* this is Danielson-Lanczos formula */
                j=i+mmax;
                tempr=wr*data[j]-wi*data[j+1];
                tempi=wr*data[j+1]+wi*data[j];
                data[j]=data[i]-tempr;
                data[j+1]=data[i+1]-tempi;
                data[i] += tempr;
                data[i+1] += tempi;
            }
            /* trigonometric recurrence */
            wr=(wtemp=wr)*wpr-wi*wpi+wr;
            wi=wi*wpr+wtemp*wpi+wi;
        }
        mmax=istep;
    }

    return;
}

/*-------------------------------------------------------------------------*
 * realft()
 *
 * from "numerical recipes in C".
 * Calculates the Fourier Transform of a set of 2*n real-valued data points.
 * Replaces this data (which is stored in the array data[1..2n]) by the
 * positive frequancy half of its complex Fourier Transform. The real-valued
 * first and last components of the complex transform are returned as elements
 * data[1] and data[2] respectively. n must be a power of 2. This routine
 * also calculates the inverse transform of a complex data array if it is the
 * tranform of real data. (Results in this case must be multiplied by 1/n.)
 *--------------------------------------------------------------------------*/

void realft(
    float *data,               /* i/o: data array   .......... */
    short n,                   /* i  : length of data array    */
    short isign                /* i  : sign +1 or -1           */
)
{
    short i,i1,i2,i3,i4,n2p3;
    float c1=0.5,c2,h1r,h1i,h2r,h2i;
    float wr,wi,wpr,wpi,wtemp,theta;

    theta=(float) (EVS_PI/(float) n);
    if (isign == 1)
    {
        /* the forward transorm here */
        c2 = -0.5;
        four1(data,n,1);
    }
    else
    {
        /* otherwise set up for the inverse transform */
        c2=0.5;
        theta = -theta;
    }
    wtemp=(float) sin(0.5f*theta);
    wpr = -2.0f*wtemp*wtemp;
    wpi=(float) sin(theta);
    wr=1.0f+wpr;
    wi=wpi;
    n2p3=2*n+3;
    /* case i=1 done separately below */
    for (i=2; i<=n/2; i++)
    {
        i4=1+(i3=n2p3-(i2=1+(i1=i+i-1)));
        /* the two separate transforms are separated out of data */
        h1r=c1*(data[i1]+data[i3]);
        h1i=c1*(data[i2]-data[i4]);
        h2r = -c2*(data[i2]+data[i4]);
        h2i=c2*(data[i1]-data[i3]);
        /* here they are recombined to form the true transform
        of the original real data */
        data[i1]=h1r+wr*h2r-wi*h2i;
        data[i2]=h1i+wr*h2i+wi*h2r;
        data[i3]=h1r-wr*h2r+wi*h2i;
        data[i4] = -h1i+wr*h2i+wi*h2r;
        /* the recurrence */
        wr=(wtemp=wr)*wpr-wi*wpi+wr;
        wi=wi*wpr+wtemp*wpi+wi;
    }
    if (isign == 1)
    {
        /* squeeze the first and the last data together to get them
        all within the original data */
        data[1] = (h1r=data[1])+data[2];
        data[2] = h1r-data[2];
    }
    else
    {
        /* this is the inverse transform for the case isign=-1 */
        data[1]=c1*((h1r=data[1])+data[2]);
        data[2]=c1*(h1r-data[2]);
        four1(data,n,-1);
    }

    return;
}