/*--------------------------------------------------------------------------
 21kcode/fft/rad4/twidrad4.c
 C program to generate the radix-4 twiddle factor table for the
 fast FFT-program RAD4.ASM for the ADSP21020.

 Based on MATLAB-File by Karl Schwarz, Erlangen-Nuernberg Universitaet

 Author  12-MAR-1991	  Ronnin Yee, Analog Devices
 Revised 25-MAR-1991	  Steven Cox, Analog Devices
---------------------------------------------------------------------------*/

#include <stdio.h>
#include <math.h>
#define	MAX_LEN	16384

int b[MAX_LEN];
main() {
    double tc;
    double ts;
    double  k1, k2, k3;		    /* coeff for table calculation      */
    double  pi;
    char    cosfn[20];  /* strings for file names */
    char    sinfn[20];
    int	    length,		    /* length of FFT */
	    iter,		    /* length/4	     */
	    i, j, k;
    FILE *s_file, *c_file;


/* initialize pi */
    pi=4.0*atan(1.0);

/* User interface */

    printf("%c%c%c%c",27,91,50,74);	/* clears ibm screen */
    printf("%c%c%c",27,91,72);		/* homes ibm cursor  */

printf("\n	 		    FFTR4TBL \n");
printf("	 Cosine and Sine Table generator for FFTRAD4.ASM \n");
printf("\nThis program generates the cosine and sine tables for the fast\n");
printf("FFT program FFTRAD4.ASM on the ADSP-21020.  The length of the FFTs\n");
printf("can be any power of four equal or greater than 64.  This program is\n");
printf("configured to generate tables for FFTs as long as %d.\n\n",MAX_LEN);

    printf("\n Enter the FFT length (max %d): ",MAX_LEN);
    scanf(" %d",&length);

    printf("\n Name of Cosine table to be created: ");
    scanf(" %s",cosfn);
    c_file=fopen(cosfn,"w");

    printf(" Name of Sine table to be created: ");
    scanf(" %s",sinfn);
    s_file=fopen(sinfn,"w");
    
/* Start Calculations	*/
    if (length > 1024)
	printf("\n\n	    Thinking hard . . .");
    else
	printf("\n\n	    Thinking . . .");

    iter = length/4;

    /* generate array for bit reversed addressing */
    for (i=1,j=1;i<=iter;i++) {
	b[i-1] = j-1;
	k = iter/2;
	while (k<j && k!=0) {
	    j =j-k;
	    k = k/2;
	}
	j += k;
    }

    /* calculate tables */

    k1 = (2*pi/(double)length);
    k2 = (2*pi/(double)length)*2;
    k3 = (2*pi/(double)length)*3;

    for (i=0;i<iter;i++) {
      tc = cos(b[i]*k1);
      ts = sin(b[i]*k1);
      fprintf(c_file,"%22.14e\n",tc);
      fprintf(s_file,"%22.14e\n",ts);

      tc = cos(b[i]*k3);
      ts = sin(b[i]*k3);
      fprintf(c_file,"%22.14e\n",tc);
      fprintf(s_file,"%22.14e\n",ts);

      tc = cos(b[i]*k2);
      ts = sin(b[i]*k2);
      fprintf(c_file,"%22.14e\n",tc);
      fprintf(s_file,"%22.14e\n",ts);
    }
    
   fclose(c_file);
   fclose(s_file);
   printf("\n    Done!\n\n");
}