/*********************************************************************\
*                                                                     *
*                             F F T 2 R 2 . C                         *
*                             ===============                         *
*                                                                     *
*                RUTINA PARA CALCULO DE LA TRANSFORMADA DE            *
*                       FOURIER EN DOS DIMENSIONES                    *
*                          ===================                        *
*                                                                     *
*                Sirve para un n£mero de puntos que sea               *
*                potencia de 2.                                       *
*                                                                     *
*                Utiliza matrices de una s¢la dimensi¢n.              *
*                                                                     *
*                                                                     *
\*********************************************************************/

# include <stdio.h>
# include <stdlib.h>
# include <math.h>

void fft2r2 (int nfil, int ncol, int isign, float x[], float y[]);
void fftr2 (int l, int n, int s, float *x, float *y, int *inte, float *seno, float *cose);


void fft2r2 (nfil, ncol, isign, x, y)

int    nfil, ncol, isign;
float  x[], y[];
{



	float    *a, *b;
	float    *seno, *cose;
	float    rnfil, rncol, rntot, var, fd, app;
	int      nf, nc, mf, mc, ld;
	int      inte[16];
	long int   i, j, ii, jj, k, ik;

	nf = nfil;   if(nfil < ncol) nf = ncol;
	a=(float *)calloc(nfil*ncol,sizeof(float));
	b=(float *)calloc(nfil*ncol,sizeof(float));
	seno=(float *)calloc(nfil*ncol,sizeof(float));
	cose=(float *)calloc(nfil*ncol,sizeof(float));

/*      ----   c lculo de la Transformada de Fourier   ----          */

	nf = nfil / 2;
	nc = ncol / 2;
	rnfil = 1./sqrt ((float)nfil);
	rncol = 1./sqrt ((float)ncol);
	rntot = rnfil * rncol;
	mf = log10 ((float)nfil) / .301030 + .1;
	mc = log10 ((float)ncol) / .301030 + .1;

/*      ----            Transformada por filas          ----          */

    /***************************************************\
    *   Esta parte se ha copiado de la fft.c para       *
    *   optimizar el tiempo de c lculo                  *
    \***************************************************/

	fd = (float)isign * 6.283185 / (float)ncol;
	ld = ncol;
	for (i = 0; i < mc; i++)        {
		ld = ld / 2;
		inte[i] = ld;
	}
	for (i = 0; i < ncol; i++)        {
	   app = (float)i * fd;
	   cose[i] = cos (app);
	   seno[i] = sin (app);
	}

    /****************************************************/

	for (ii = 0; ii < nfil; ii++)   {
		k = ii * (long)ncol;
		for (i = 0; i < nc; i++)  {
			ik = k + i;
			j  = i + nc;
			jj = k + j;
			a[i] = x[jj];
			b[i] = y[jj];
			a[j] = x[ik];
			b[j] = y[ik];
		}

		fftr2 (ncol, mc, isign, a, b, inte, seno, cose);

		for (i = 0; i < nc; i++)  {
			ik = k + i;
			j  = i + nc;
			jj = k + j;
			x[jj] = a[i] ;
			y[jj] = b[i] ;
			x[ik] = a[j] ;
			y[ik] = b[j] ;
		}
	}


/*      ----            Transformada por columnas       ----          */


    /***************************************************\
    *   Esta parte se ha copiado de la fft.c para       *
    *   optimizar el tiempo de c lculo                  *
    \***************************************************/

	if ( nfil != ncol ) {
	  fd = (float)isign * 6.283185 / (float)nfil;
	  ld = nfil;
	  for (i = 0; i < mf; i++)        {
		  ld = ld / 2;
		  inte[i] = ld;
	  }
	  for (i = 0; i < nfil; i++)        {
	      app = (float)i * fd;
	      cose[i] = cos (app);
	      seno[i] = sin (app);
	  }
	}

    /*************************************************/
	k = nf * (long) ncol;
	for (jj = 0; jj < ncol; jj++)   {
		for (i = 0; i < nf; i++)  {
			j = i + nf;
			ii = i * ncol + jj;
			ik = ii + k;
			a[i] = x[ik];
			b[i] = y[ik];
			a[j] = x[ii];
			b[j] = y[ii];
		}

		fftr2 (nfil, mf, isign, a, b, inte, seno, cose);

		for (i = 0; i < nf; i++)  {
			j = i + nf;
			ii = i * ncol + jj;
			ik = ii + k;
			x[ik] = a[i] * rntot;
			y[ik] = b[i] * rntot;
			x[ii] = a[j] * rntot;
			y[ii] = b[j] * rntot;
		}
	}
	return;
}



 /********************************************************************\
 *                                                                    *
 *                        F F T   . C                                 *
 *                        ===========                                 *
 *                                                                    *
 *      Funci¢n para calcular la transformada de Fourier.             *
 *      Trabaja con un n£mero de puntos que sea potencia de 2.        *
 *      Par metros:                                                   *
 *              l       n£mero de puntos                              *
 *              n       exponente de la potencia de 2                 *
 *              s         1 para transformada directa                 *
 *                      - 1 para transformada inversa                 *
 *              x       vector de datos ( parte real )                *
 *              y       vector de datos ( parte imaginaria )          *
 *                                                                    *
 *                                                                    *
 \********************************************************************/



void fftr2 (l, n, s, x, y, inte, seno, cose)
	int     l, n, s, *inte;
	float    *x, *y;
	float    *seno, *cose;
{
	int     ld, lt, lc, lci, i, m;
	int     nt, nd, nc, ib, j, k, id, nci;
	float   fd, fc, qu, qd,  hu, hd;

	nd = 1;
	for (m = 1; m <= n; m++)        {
		nt = nd;
		nd = nd + nd;
		lt = l / nt;
		lc = lt / 2;
		nc = 0;
		for (ib = 0; ib < nt; ib++)    {
			lci = lt * ib;
			for (i = 0; i < lc; i++)      {
				j = i + lci;
				k = j + lc;
				qu = x[k] * cose[nc] - y[k] * seno[nc];
				qd = x[k] * seno[nc] + y[k] * cose[nc];
				x[k] = x[j] - qu;
				y[k] = y[j] - qd;
				x[j] = x[j] + qu;
				y[j] = y[j] + qd;
			}
			for (i = 1; i < n; i++) {
				id = i;
				if (inte[i] > nc)        goto r40;
				nc = nc - inte[i];
			}
r40:                    nc = nc + inte[id];

		}
	}
	nc = 0;
	for (k = 0; k < l; k++) {
		if (nc > k)     {
			hu = x[nc];
			hd = y[nc];
			x[nc] = x[k];
			y[nc] = y[k];
			x[k] = hu;
			y[k] = hd;
		}
		for (i = 0; i < n; i++)     {
			id = i;
			if (inte[i] > nc)    goto r80;
			nc = nc - inte[i];
		}
r80:            nc = nc + inte[id];
	}
	return;
}