/*--------------------------------------------------------------------------
 21kcode/fft/rad4/irf4tbl.c
 C program to generate the radix-4 twiddle factor table for the
 fast FFT-program IRFFT4.ASM for the ADSP-21020.

 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
         15-APR-1991	  Ronnin Yee, Analog Devices
	 12-JUN-1991	  Ronnnin Yee, Analog Devices
---------------------------------------------------------------------------*/

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

main() {
    double tc,c;
    double ts,s;
    double  k1, k2, k3;		    /* coeff for table calculation      */
    double  pi, freq;
    char    cosfn[20];  /* strings for file names */
    char    sinfn[20];
    int	    *b;			    /* pntr to bitrev array */
    int	    length,		    /* length of FFT */
	    iter,		    /* length/4	     */
	    i, j, k, n;
    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	 		    IRF4TBL \n");
printf("	 Cosine and Sine Table generator for IRFFT4.ASM \n");
printf("\nThis program generates the cosine and sine tables for the inverse\n");
printf("RFFT routine IRFFT4.ASM on the ADSP-21020.  The length of the IRFFTs\n");
printf("can be twice any power of four and must be equal or greater than 128.\n");
printf("Note that the twiddle factor tables are identical to the ones used by\n");
printf("RFFT4.ASM.  \n");

    printf("\n Enter the FFT length: ");
    scanf(" %d",&length);
    length=length/2;

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

    printf(" Name of Sine Twiddle Factor table to be created: ");
    scanf(" %s",sinfn);
    s_file=fopen(sinfn,"w");

/* allocate the storage for bitrev array */
    b = (int *) calloc(length,sizeof(int));
    if(!b) {
        printf("\nUnable to allocate memory for bitrev. array.\n");
        exit(1);
    }
    
/* 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);

/* Make (1/N)*cos(pi*n/N) and (1/N)*sin(pi*n/N) tables */

  n=length;
  iter=n/2-1;

  printf("\n\n  Enter the (1/N)*cosine Conversion Coeff table filename: ");
  scanf(" %s",cosfn);
  c_file=fopen(cosfn,"w");

  printf("  Enter the (1/N)*sine Conversion Coeff table filename: ");
  scanf(" %s",sinfn);
  s_file=fopen(sinfn,"w");

  freq=pi/(double)n;
  k1=1.0/(2.0*n);
  for (i=0; i <= iter; i++)
  {
    s=k1*sin((double)i * freq);
    c=k1*cos((double)i * freq);
    fprintf(s_file,"%22.14e\n",s);
    fprintf(c_file,"%22.14e\n",c);
    printf("%d : %22.14e : %22.14e\r",i,s,c);
  }
  fclose(s_file);
  fclose(c_file);

   printf("\n\n    Done!\n\n");
}