/*--------------------------------------------------------------------------
 21kcode/fir/ffir/ffircoef.c
 C program to generate complex coefficients for FFIR.ASM from ordinary 
 direct-form FIR filter taps.  The coefficients are = FFT(taps+zeros).

 Author  10-JUN-1991	  Ronnin Yee, Analog Devices
---------------------------------------------------------------------------*/


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

/* COMPLEX STRUCTURE */
typedef struct {
    double real, imag;
} COMPLEX;

void fft();

#define LINELEN    256

main () 
{
    FILE    *infile,*outfile;
    char    line[LINELEN],filename[25];
    int	    len, points, stages, i;
    double  tmp_double;
    COMPLEX *data,*tmp_data;


  printf("%c%c%c%c",27,91,50,74); /*clear screen*/
  printf("%c%c%c",27,91,72);      /*curser home*/
  printf("______Transfer Function Coefficient Generator for FFIR.ASM_____\n");
  printf("This program converts regular FIR tap coefficients into coeff.\n");
  printf("for FFIR.ASM.  After querying for the size of the FFTs to be \n");
  printf("done in FFIR.ASM, this program reads the FIR coefficients from\n");
  printf("a file, adds zeros to pad them to the FFT length, performs a\n");
  printf("FFT on them and saves the real and imaginary results in files.\n");
  printf("\nVersion 1.0, 9-JUN-91, Ronnin Yee, Analog Devices\n");
  printf("_______________________________________________________________\n");

    printf("\nEnter the name of the file of FIR coeffiecients: ");
    scanf(" %s",filename);
    if((infile = fopen(filename, "r")) == NULL) {
	fprintf(stderr, "Error opening %s\n", filename);
	exit(1);
    }

/* scan incoming file */
    len=0;
    while (fgets(line, LINELEN, infile) != NULL) {
	len++;
	if (sscanf(line, "%lf", &tmp_double)!=1) {
	   fprintf(stderr, 
		"Error: floating point number not found in %s\n",filename);
	   fprintf(stderr, "       line %d:%s\n",len,line);
	   exit(1);
	}
    }

    rewind(infile);

/* get FFT size */
    printf("\n    %d FIR filter coefficients have been read in.",len);
  printf("\n    The suggested FFT size is %d points for %d taps.\n",3*len,len);
    printf("\nEnter the size of the FFTs to be done: ");
    scanf(" %d",&points);
    if (points <= len){
	fprintf(stderr,"Error: Number of Taps must be less than the FFT size");
	exit(1);
    }

/* allocate the storage for data */
    data = (COMPLEX *) calloc(points,sizeof(COMPLEX));
    if(!data) {
        printf("\nUnable to allocate complex data array\n");
        exit(1);
    }

/* load up data array */
/*	note that all unwritten values of real and all imag are zero	*/
    tmp_data=data;
    for (i=0;i<len;i++) {
	fgets(line, LINELEN, infile);
	sscanf(line, "%lf", &tmp_double);
	tmp_data->real = tmp_double;
	tmp_data++;
    }
    fclose(infile);

/* do the fft */
    stages=(int)(log((double)points)/log((double)2.0));
    fft(data,stages);

/* write out the coefficients */
    printf("\nEnter the output filename for real coefficients: ");
    scanf(" %s",filename);
    outfile=fopen(filename,"w");

    tmp_data=data;
    for (i=0;i<points;i++) {
	fprintf(outfile,"%22.14e\n",tmp_data->real);
	tmp_data++;
    }
    fclose(outfile);

    printf("\nEnter the output filename for imaginary coefficients: ");
    scanf(" %s",filename);
    outfile=fopen(filename,"w");

    tmp_data=data;
    for (i=0;i<points;i++) {
	fprintf(outfile,"%22.14e\n",tmp_data->imag);
	tmp_data++;
    }
    fclose(outfile);

    free(data);
    printf("\n");
}


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

fft - In-place radix 2 decimation in time FFT

Requires pointer to complex array, x and power of 2 size of FFT, m
(size of FFT = 2**m).  Places FFT output on top of input COMPLEX array.

void fft(COMPLEX *x, int m)

This routine has been modified from the fft routine in
"C Language Algorithms for Digital Signal Processing" by Embree & Kimble

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

void fft(x,m)
    COMPLEX *x;
    int m;
{
    static COMPLEX *w;           /* used to store the w complex array */
    static int mstore = 0;       /* stores m for future reference */
    static int n = 1;            /* length of fft stored for future */

    COMPLEX u,temp,tm;
    COMPLEX *xi,*xip,*xj,*wptr;

    int i,j,k,l,le,windex;

    double arg;

    if(m != mstore) {

/* free previously allocated storage and set new m */

        if(mstore != 0) free(w);
        mstore = m;
        if(m == 0) return;       /* if m=0 then done */

/* n = 2**m = fft length */

        n = 1 << m;
        le = n/2;

/* allocate the storage for w */

        w = (COMPLEX *) calloc(le-1,sizeof(COMPLEX));
        if(!w) {
            printf("\nUnable to allocate complex W array\n");
            exit(1);
        }

/* calculate the w values (double precision) */

        arg = 4.0*atan(1.0)/ (double)le;         /* PI/le calculation */
        xj = w;
        for (j = 1 ; j < le ; j++) {
            xj->real = cos((double)j * arg);
            xj->imag = -sin((double)j * arg);
            xj++;
        }
    }

/* start fft */

    le = n;
    windex = 1;
    for (l = 0 ; l < m ; l++) {
        le = le/2;

/* first iteration with no multiplies */

        for(i = 0 ; i < n ; i = i + 2*le) {
            xi = x + i;
            xip = xi + le;
            temp.real = xi->real + xip->real;
            temp.imag = xi->imag + xip->imag;
            xip->real = xi->real - xip->real;
            xip->imag = xi->imag - xip->imag;
            *xi = temp;
        }

/* remaining iterations use stored w */

        wptr = w + windex - 1;
        for (j = 1 ; j < le ; j++) {
            u = *wptr;
            for (i = j ; i < n ; i = i + 2*le) {
                xi = x + i;
                xip = xi + le;
                temp.real = xi->real + xip->real;
                temp.imag = xi->imag + xip->imag;
                tm.real = xi->real - xip->real;
                tm.imag = xi->imag - xip->imag;             
                xip->real = tm.real*u.real - tm.imag*u.imag;
                xip->imag = tm.real*u.imag + tm.imag*u.real;
                *xi = temp;
            }
            wptr = wptr + windex;
        }
        windex = 2*windex;
    }            

/* rearrange data by bit reversing */

    j = 0;
    for (i = 1 ; i < (n-1) ; i++) {
        k = n/2;
        while(k <= j) {
            j = j - k;
            k = k/2;
        }
        j = j + k;
        if (i < j) {
            xi = x + i;
            xj = x + j;
            temp = *xj;
            *xj = *xi;
            *xi = temp;
        }
    }
}