Files
rtphone/src/libs/libevs/lib_com/ifft_rel.cpp
2020-06-15 15:41:54 +03:00

256 lines
6.4 KiB
C++

/*====================================================================================
EVS Codec 3GPP TS26.443 Nov 13, 2018. Version 12.11.0 / 13.7.0 / 14.3.0 / 15.1.0
====================================================================================*/
#include "options.h"
#include "prot.h"
#include "rom_com.h"
/*-----------------------------------------------------------------*
* Local constants
*-----------------------------------------------------------------*/
#define N_MAX_FFT 1024
#define INV_SQR2 0.70710676908493f
/*---------------------------------------------------------------------*
* ifft_rel()
*
* Calculate the inverse FFT of a real signal
*
* Based on the FORTRAN code from the article "Real-valued Fast Fourier Transform Algorithms"
* by Sorensen, ... in IEEE Trans. on ASSP, Vol. ASSP-35, No. June 6th 1987.
*
* Input: the io[] signal containing the spectrum in the following order :
*
* Re[0], Re[1], .. Re[n/2], Im[n/2-1], .. Im[1]
*---------------------------------------------------------------------*/
void ifft_rel(
float io[], /* i/o: input/output vector */
const short n, /* i : vector length */
const short m /* i : log2 of vector length */
)
{
short i, j, k;
short step;
short n2, n4, n8, i0;
short is, id;
float *x,*xi0, *xi1, *xi2, *xi3, *xi4, *xup1, *xdn6, *xup3, *xdn8;
float xt;
float r1;
float t1, t2, t3, t4, t5;
float cc1, cc3, ss1, ss3;
const float *s, *s3, *c, *c3;
const short *idx;
float temp[512];
float n_inv;
n_inv = 1.0f/n;
/*-----------------------------------------------------------------*
* IFFT
*-----------------------------------------------------------------*/
x = &io[-1];
n2 = 2*n;
for (k=1; k<m; k++)
{
is = 0;
id = n2;
n2 = n2 >> 1;
n4 = n2 >> 2;
n8 = n4 >> 1;
while (is < n-1)
{
xi1 = x + is + 1;
xi2 = xi1 + n4;
xi3 = xi2 + n4;
xi4 = xi3 + n4;
for (i=is; i<n; i+= id)
{
t1 = *xi1 - *xi3;
*xi1 += *xi3;
*xi2 = 2.0f**xi2;
*xi3 = t1- 2.0f**xi4;
*xi4 = t1 + 2.0f**xi4;
if (n4 != 1)
{
t1 = (*(xi2+n8) - *(xi1+n8))*INV_SQR2;
t2 = (*(xi4+n8) + *(xi3+n8))*INV_SQR2;
*(xi1+n8) += *(xi2+n8);
*(xi2+n8) = *(xi4+n8) - *(xi3+n8);
*(xi3+n8) = (float)(2.0f * (-t2-t1));
*(xi4+n8) = (float)(2.0f * (-t2+t1));
}
xi1 += id;
xi2 += id;
xi3 += id;
xi4 += id;
}
is = 2*id - n2;
id = 4*id;
}
step = N_MAX_FFT/n2;
s = sincos_t_ext + step;
c = s + N_MAX_FFT/4;
s3 = sincos_t_ext + 3*step;
c3 = s3 + N_MAX_FFT/4;
for (j=2; j<=n8; j++)
{
cc1 = *c ;
ss1 = *s;
cc3 = *c3;
ss3 = *s3;
is = 0;
id = 2 * n2;
c += step;
s += step;
c3 += 3*step;
s3 += 3*step;
while (is < n-1)
{
xup1 = x + j + is;
xup3 = xup1 + 2*n4;
xdn6 = xup3 - 2*j +2;
xdn8 = xdn6 + 2*n4;
for (i=is; i<n; i+=id)
{
t1 = *xup1 - *xdn6;
*xup1 = *xup1 + *xdn6;
xup1 += n4;
xdn6 -= n4;
t2 = *xdn6 - *xup1;
*xdn6 = *xup1 + *xdn6;
xdn6 += n4;
t3 = *xdn8 + *xup3;
*xdn6 = *xdn8 - *xup3;
xup3 += n4;
xdn8 -= n4;
t4 = *xup3 + *xdn8;
*xup1= *xup3 - *xdn8;
t5 = t1 - t4;
t1 = t1 + t4;
t4 = t2 - t3;
t2 = t2 + t3;
*xup3 = t1*cc3 - t2*ss3;
xup3 -= n4;
*xup3 = t5*cc1 + t4*ss1;
*xdn8 = -t4*cc1 + t5*ss1;
xdn8 += n4;
*xdn8 = t2*cc3 + t1*ss3;
xup1 -= n4;
xup1 += id;
xup3 += id;
xdn6 += id;
xdn8 += id;
}
is = 2*id - n2;
id = 4*id;
}
}
}
/*-----------------------------------------------------------------*
* Length two butterflies
*-----------------------------------------------------------------*/
is = 1;
id = 4;
while (is < n)
{
xi0 = x + is ;
xi1 = xi0 + 1;
for (i0=is; i0<=n; i0+=id)
{
r1 = *xi0;
*xi0= r1 + *xi1;
*xi1 = r1 - *xi1;
xi0 += id;
xi1 += id;
}
is = 2*id - 1;
id = 4*id;
}
/*-----------------------------------------------------------------*
* Digit reverse counter
*-----------------------------------------------------------------*/
idx = fft256_read_indexes;
xi0 = temp-1;
if (n == 128)
{
for (i=0; i<n; i++)
{
j = *idx++;
temp[i] = x[1+(j>>1)];
}
}
else if (n == 256)
{
for (i=0; i<n; i++)
{
j = *idx++;
temp[i] = x[1+j];
}
}
else if (n == 512)
{
for (i=0; i<256; i++)
{
j = *idx++;
temp[i] = x[1+2*j];
temp[i+256] = x[2+2*j];
}
}
else
{
xi0 = x;
j = 1;
for (i=1; i<n; i++)
{
if (i < j)
{
xt = x[j];
x[j] = x[i];
x[i] = xt;
}
k = n >> 1;
while (k < j)
{
j = j - k;
k = k >> 1;
}
j = j + k;
}
}
/*-----------------------------------------------------------------*
* Normalization
*-----------------------------------------------------------------*/
for (i=1; i<=n; i++)
{
x[i] = xi0[i] * n_inv;
}
return;
}