#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <graphics.h>
#include "freq.h"

int fftlen=1024;        // Number of points for FFT
long SampleRate=44100;  // A/D sampling rate
int logfreq=0;          // Flag set to 1 for log-based frequency scale
int logamp=0;           // Flag set to 1 for log-based amplitude scale
int windfunc=0;         // Flag set to selected window function
int ys=10;              // Flag set for max of y-axis (default=10 = 1.0)
int logb=6;             // Flag set for base of log scale (default=6 = -60db)
int logs=0;             // Flag for max of log scale (default=0 = 0db)
int gain=0;             // Amount of db/octave gain
int gain3db=0;          // Flag indicating a 3db/octave scale factor gain
int deriv=0;            // Flag for doing differencing for 6db/octave gain
long ref_freq=1000;     // Reference frequency for n db/octave gains

char *windfile=NULL;    // Pointer to filename for saving window data

extern int *buffer[];   // Buffers for gathering data
extern int *fftdata;    // Array for FFT data
extern int *wind;       // Array storing windowing function

extern int flag[];            // Array of flags indicating fullness of buffers
extern int x[];               // Array of bin #'s displayed
extern int lasty[];           // Last y position for FFT bins
extern unsigned int yscale[]; // Scaling factors
extern int ybase[];           // Scaling offset for log calculations
extern int shift;             // Number of bits for gain shift


double alpha;           // Gaussian window parameter

struct rgb background = { 0,0,20 },
           warn       = { 20,0,0 },
           graph      = { 30,35,60 },
           tick       = { 40,40,40 },
           label      = { 50,20,45 },
           border     = { 40,40,40 },
           text       = { 55,55,25 },
           darkhl     = { 63,63,63 },
           lighthl    = { 55,55,55 };

/*
 *  Parse the ini file, if it exists
 */
void parse_ini_file(void)
{
   FILE *fp;
   char buffer[200];

   if((fp=fopen("FREQ.INI","r"))!=NULL)
   {
      while(!feof(fp))
      {
         fgets(buffer,200,fp);

         sscanf(buffer,"Sample rate:%ld",&SampleRate);
         sscanf(buffer,"FFT length:%d",&fftlen);
         sscanf(buffer,"Window function:%d",&windfunc);
         sscanf(buffer,"Log freq scale:%d",&logfreq);
         sscanf(buffer,"Log amp scale:%d",&logamp);
         sscanf(buffer,"Base db:%d",&logb);
         sscanf(buffer,"Top db:%d",&logs);
         sscanf(buffer,"Max amp:%d",&ys);
         sscanf(buffer,"DB/octave gain:%d",&gain);
         sscanf(buffer,"Reference frequency:%d",&ref_freq);
         sscanf(buffer,"Background color:%d,%d,%d",&background.red,&background.green,&background.blue);
         sscanf(buffer,"Clipping warning color:%d,%d,%d",&warn.red,&warn.green,&warn.blue);
         sscanf(buffer,"Graph color:%d,%d,%d",&graph.red,&graph.green,&graph.blue);
         sscanf(buffer,"Tick mark color:%d,%d,%d",&tick.red,&tick.green,&tick.blue);
         sscanf(buffer,"Axis label color:%d,%d,%d",&label.red,&label.green,&label.blue);
         sscanf(buffer,"Border color:%d,%d,%d",&border.red,&border.green,&border.blue);
         sscanf(buffer,"Text color:%d,%d,%d",&text.red,&text.green,&text.blue);
         sscanf(buffer,"Cursor upper color:%d,%d,%d",&darkhl.red,&darkhl.green,&darkhl.blue);
         sscanf(buffer,"Cursor lower color:%d,%d,%d",&lighthl.red,&lighthl.green,&lighthl.blue);
      }
      fclose(fp);
   }
}


/*
 * Parse the command line
 */
void parse_command(int argc,char *argv[])
{
   int i=1;
   while(i<argc)
   {
      if(argv[i][0]=='-')
      {
         switch(argv[i][1])
         {
            case 'S':
            case 's':
                 SampleRate=atol(&argv[i][2]);
                 break;
            case 'F':
            case 'f':
                 if((fftlen=atoi(&argv[i][2])) > MAX_LEN)
                 {
                    printf("Maximum of %d points allowed!\n",MAX_LEN);
                    puts("Hit any key to continue with the maximum...");
                    fftlen=MAX_LEN;
                    getch();
                 }
                 break;
            case 'M':
            case 'm':
                 ys=atoi(&argv[i][2]);
                 break;
            case 'T':
            case 't':
                 logs=atoi(&argv[i][2]);
                 break;
            case 'B':
            case 'b':
                 logb=atoi(&argv[i][2]);
                 break;
            case 'G':
            case 'g':
                 gain=atoi(&argv[i][2]);
                 break;
            case 'R':
            case 'r':
                 ref_freq=atol(&argv[i][2]);
                 break;
            case 'L':
            case 'l':
                 if((argv[i][2]=='A')||(argv[i][2]=='a'))
                    logamp=1;
                 else if((argv[i][2]=='F')||(argv[i][2]=='f'))
                    logfreq=1;
                 else
                 {
                    printf("Ignoring unrecognized switch: %s\n",argv[i]);
                    puts("Hit any key to continue...");
                    getch();
                 }
                 break;
            case 'W':
            case 'w':
                 if(argv[i][2]>='A')
                 {
                    windfile=&argv[i][2];
                 }
                 else
                 {
                    windfunc=argv[i][2]-'0';
                    if(windfunc==4)
                    {
                       alpha=0;
                       if(argv[i][3]==',')
                          alpha=atof(&argv[i][4]);
                       if(alpha<=0) alpha=1;
                    }
                 }

                 break;
            case '?':
            case 'H':
            case 'h':
                 puts("-Snumber sets the sampling rate.");
                 puts("-Fnumber sets the length of the FFT.");
                 puts("-Mnumber sets the scale maximum.");
                 puts("-Bnumber sets the logarithmic base level (in tens of dB).");
                 puts("-Tnumber sets the logarithmic top level (in tens of dB).");
                 puts("-Gnumber sets the dB/octave gain factor.");
                 puts("-Rnumber sets the reference frequency for dB/octave gains.");
                 puts("-LA sets a logarithmic scale for the amplitude axis.");
                 puts("-LF sets a logarithmic scale for the frequency axis.");
                 puts("-W0 selects a Hamming window.  (offset sine) <--default");
                 puts("-W1 selects a Hanning window.  (sine)");
                 puts("-W2 selects a Blackman window. (two sines)");
                 puts("-W3[,alpha] selects a Gaussian window.");
                 puts("-W4 selects a Welch window.    (quadratic)");
                 puts("-W5 selects a Parzen window.   (triangular)");
                 puts("-W6 selects a Rectangular window.");
                 puts("-Wfilename saves the windowing function to the specified file.");
                 puts("-? or -H displays this message.");
                 exit(0);
            default:
                 printf("Ignoring unrecognized switch: %s\n",argv[i]);
                 puts("Hit any key to continue...");
                 getch();
         }
      }
      else
         printf("Ignoring unrecognized parameter: %s  (use -? for help)\n",argv[i]);
      i++;
   }

   /*
    * Set up the flags based on chosen gain
    */
   if(gain>11)
   {
      gain=12;
      gain3db=0;
      deriv=2;
   }
   else if(gain>8)
   {
      gain=9;
      gain3db=1;
      deriv=1;
   }
   else if(gain>5)
   {
      gain=6;
      gain3db=0;
      deriv=1;
   }
   else if(gain>2)
   {
      gain=3;
      gain3db=1;
      deriv=0;
   }
   else
   {
      gain=0;
      gain3db=0;
      deriv=0;
   }

   /*
    * Watch for bad choice of log scales
    */
   if(logb<=logs) logs=logb-3;

   if((SampleRate>88200) || (SampleRate<5000))
   {
      puts("Sample rate must be between 5000 and 88200");
      exit(1);
   }

}

/*
 *  Allocate memory for and initialize the buffers
 */
void setup_buffers(int length)
{
   int i;

   for(i=0;i<BUFFERS;i++)
   {
      if((buffer[i]=(int *)malloc(length*sizeof(int)))==NULL)
      {
         puts("Unable to allocate enough memory for the buffers!");
         exit(1);
      }
      flag[i]=0;
   }

   if((fftdata=(int *)malloc(length*sizeof(int)))==NULL)
   {
      puts("Unable to allocate enough memory for the buffers!");
      exit(1);
   }
   if((wind=(int *)malloc(length*sizeof(int)))==NULL)
   {
      puts("Unable to allocate enough memory for the buffers!");
      exit(1);
   }

   /*
    *  Clear out the buffer of last Y display positions
    */
   for(i=0;i<(WINDOW_RIGHT-WINDOW_LEFT+1);i++)
      lasty[i]=WINDOW_BOTTOM;
}

void compute_window_function(void)
{
   int i;
   double val;

   /*
    *  Calculate FFT Windowing function
    */
   for(i=0;i<fftlen;i++)
   {
      double val;

      switch(windfunc)
      {
         case 1:    // Hanning
            val = 0.5 - 0.5*cos(2*M_PI*i/fftlen);
            break;

         case 2:    // Blackman
            val = 0.42-0.5*cos(2*M_PI*i/fftlen)+0.08*cos(4*M_PI*i/fftlen);
            break;

         case 3:    // Gaussian
            val = exp(-alpha / ((double)fftlen*(double)fftlen)
                             * ((double)2*i-fftlen)*((double)2*i-fftlen));
            break;

         case 4:    // Welch
            val = 1 -  ((double)(2*i-fftlen)/(double)(fftlen+1))
                      *((double)(2*i-fftlen)/(double)(fftlen+1));
            break;

         case 5:    // Parzen
            val = 1 - fabs((double)(2*i-fftlen)/(double)(fftlen+1));
            break;

         case 6:    // Rectangular
            val = 1;
            break;

         default:   // Hamming
            val = 0.54-0.46*cos(2*M_PI*i/fftlen);
            break;
      }
      wind[i]=floor(val*32767+0.5);
   }

   /*
    *  If output requested, save the window function. (for the curious)
    */
   if(windfile)
   {
      FILE *fp;
      if((fp=fopen(windfile,"w"))==NULL)
      {
         printf("Error opening window output file: %s\n",windfile);
         getch();
      }
      else
      {
         for(i=0;i<fftlen;i++)
            fprintf(fp,"%d\n",wind[i]);
         fclose(fp);
      }
      windfile=NULL;
   }
}

/*
 *  Set up X axis scales
 */
void setup_xscale(void)
{
   int i;
   /*
    *  Initialize graph x scale (linear or logarithmic).
    *  This array points to the bin to be plotted on a given line.
    */
   for(i=0;i<(WINDOW_RIGHT-WINDOW_LEFT+1);i++)
   {
      if(logfreq)
         x[i]=floor(exp((double)i/(double)(WINDOW_RIGHT-WINDOW_LEFT+1)*log(fftlen/2))+0.5);
      else
         x[i]=floor((double)i*(double)fftlen/(double)(2*(WINDOW_RIGHT-WINDOW_LEFT+1))+0.5);
   }
   /*
    *  If lines are repeated on the screen, flag this so that we don't
    *  have to recompute the y values.
    */
   for(i=(WINDOW_RIGHT-WINDOW_LEFT);i>0;i--)
   {
      if(x[i]==fftlen/2)  /* we don't compute this bin, so don't display it */
         x[i]=-1;
      if(x[i]==x[i-1])
         x[i]=-1;
   }
}

/*
 *  Set up logarithmic amplitude scale factors and offsets.
 */
void setup_logscales(void)
{
   int i;
   double scale,base,convert;

   /*
    *  Compute the (logarithmic) y scale factor and offset.
    *  This may include a 3dB/octave gain.
    *
    *  Conversion factor from db/10 to dPhils (the computed "unit")
    *  where 1024 dPhils equals 2.0
    */
   convert=1024.0*log(10)/log(2);

   scale=(WINDOW_BOTTOM-WINDOW_TOP)/(logb-logs)/convert;
   if(deriv==0)
      base=(log10(32768.0)*2.0-logb)*convert;
   else if(deriv==1)
      base=(log10(32768.0)*2.0-log10(SampleRate/(2*M_PI*ref_freq))*2.0-logb)*convert;
   else
      base=(log10(32768.0)*2.0-log10(SampleRate/(2*M_PI*ref_freq))*4.0-logb)*convert;


   shift=0;
   /*
    *  Make maximum use of available bits
    */
   while(scale<32767.75)
   {
      scale*=2;
      shift++;
   }
   scale=floor(scale+0.5);

   if(gain3db)
   {
      for(i=0;i<(WINDOW_RIGHT-WINDOW_LEFT+1);i++)
      {
         if(x[i]==-1)
         {
            yscale[i]=0;
            ybase[i]=0;
         }
         else
         {
            yscale[i]=scale;
            if(x[i]!=0)
               ybase[i]=floor(base-log10((double)x[i]*SampleRate/fftlen/ref_freq)*convert+0.5);
            else
            {
               yscale[i]=0;
               ybase[i]=base;
            }
         }
      }
   }
   else
   {
      for(i=0;i<(WINDOW_RIGHT-WINDOW_LEFT+1);i++)
      {
         if(x[i]==-1)
         {
            yscale[i]=0;
            ybase[i]=0;
         }
         else
         {
            yscale[i]=scale;
            ybase[i]=floor(0.5+base);
         }
      }
   }
}

/*
 *  Set up linear amplitude scale factors
 */
void setup_linscales(void)
{
   int i;
   double scale;

   /*
    *  Compute the (linear) y scale factor.
    *  This may include a 3dB/octave gain.
    */
   scale=(WINDOW_BOTTOM-WINDOW_TOP)/(ys*3276.8*sqrt(ref_freq));
   shift=0;

   if(deriv==1)
   {
      scale*=(double)SampleRate/(2*M_PI*(double)ref_freq);
   }
   else if(deriv==2)
   {
      scale*=(double)SampleRate*(double)SampleRate
             /(4*M_PI*M_PI*(double)ref_freq*(double)ref_freq);
   }

   /*
    *  Make maximum use of available bits
    */
   if(gain3db)
   {
      while((scale*sqrt(SampleRate/2))<32767.75)
      {
         scale*=2;
         shift++;
      }
      for(i=0;i<(WINDOW_RIGHT-WINDOW_LEFT+1);i++)
      {
         if(x[i]==-1)
            yscale[i]=0;
         else
            yscale[i]=floor(0.5+scale*sqrt(x[i]*SampleRate/fftlen));
      }
   }
   else
   {
      scale=scale*sqrt(ref_freq);
      while(scale<32767.75)
      {
         scale*=2;
         shift++;
      }

      scale=floor(scale+0.5);

      for(i=0;i<(WINDOW_RIGHT-WINDOW_LEFT+1);i++)
      {
         if(x[i]==-1)
            yscale[i]=0;
         else
            yscale[i]=scale;
      }
   }
}

void frequency_scale(void)
{
   char text[20];
   /*
    *  Put up the frequency scale.
    */
   setfillstyle(SOLID_FILL,0);
   bar(WINDOW_LEFT-10,WINDOW_BOTTOM+3,WINDOW_RIGHT+10,479);

   if(logfreq)
   {
      double f,step,base;
      settextjustify(CENTER_TEXT,TOP_TEXT);
      settextstyle(0,1,0);
      step=exp(log(fftlen/2)/32.0);
      base=(double)SampleRate/(double)fftlen;
      for(f=base;f<=(double)(SampleRate*1.01)/2;f*=step)
      {
         setcolor(BORDER_COLOR);
         moveto(WINDOW_LEFT+(log(f)-log(base))/(log(SampleRate/2)-log(base))*(WINDOW_RIGHT-WINDOW_LEFT+1),WINDOW_BOTTOM+3);
         linerel(0,3);
         moverel(0,3);
         setcolor(LABEL_COLOR);
         sprintf(text,"%7.1lf",f);
         outtext(text);
      }
   }
   else
   {
      double f;
      settextjustify(CENTER_TEXT,TOP_TEXT);
      settextstyle(0,1,0);
      for(f=0;f<=(double)SampleRate/2.0;f+=(double)SampleRate/64.0)
      {
         setcolor(BORDER_COLOR);
         moveto(WINDOW_LEFT+f*(double)(WINDOW_RIGHT-WINDOW_LEFT+1)*2.0/(double)SampleRate,WINDOW_BOTTOM+3);
         linerel(0,3);
         moverel(0,3);
         setcolor(LABEL_COLOR);
         sprintf(text,"%7.1lf",f);
         outtext(text);
      }
   }
}

void amplitude_scale(void)
{
   int i;
   char text[20];

   setfillstyle(SOLID_FILL,0);
   bar(0,WINDOW_TOP-10,WINDOW_LEFT-3,WINDOW_BOTTOM+10);

   /*
    *  Put up the amplitude scale
    */
   if(logamp)
   {

      settextjustify(RIGHT_TEXT,CENTER_TEXT);
      settextstyle(0,0,0);
      for(i=0;i<(logb-logs)+1;i++)
      {
         setcolor(BORDER_COLOR);
         moveto(WINDOW_LEFT-3,WINDOW_TOP+i*(WINDOW_BOTTOM-WINDOW_TOP)/(logb-logs));
         linerel(-3,0);
         moverel(-3,0);
         setcolor(LABEL_COLOR);
         sprintf(text,"%d",(i+logs)*(-10));
         outtext(text);
      }
   }
   else
   {
      double scale=(double)(ys*3276.8)/(double)(WINDOW_BOTTOM-WINDOW_TOP);
      settextjustify(RIGHT_TEXT,CENTER_TEXT);
      settextstyle(0,0,0);
      for(i=0;i<=10;i++)
      {
         setcolor(BORDER_COLOR);
         moveto(WINDOW_LEFT-3,WINDOW_BOTTOM-i*ys*327.68/scale);
         linerel(-3,0);
         moverel(-3,0);
         setcolor(LABEL_COLOR);
         sprintf(text,"%4.2f",(float)i*ys*0.01);
         outtext(text);
      }
   }
}

char *window_name[] = { "Hamming",
                        "Hanning",
                        "Blackman",
                        "Gaussian",
                        "Welch",
                        "Parzen",
                        "Rectangular" };

void update_header(void)
{
   char ach[100];

   setfillstyle(SOLID_FILL,0);
   bar(SRX-150,SRY,SRX+150,RFY+10);

   settextstyle(0,0,0);
   settextjustify(CENTER_TEXT,TOP_TEXT);
   setcolor(TEXT_COLOR);
   sprintf(ach,"Sampling rate: %ld Hz",SampleRate);
   outtextxy(SRX,SRY,ach);
   sprintf(ach,"FFT Length: %d points",fftlen);
   outtextxy(FLX,FLY,ach);
   sprintf(ach,"Frequency resolution: %g Hz",(double)SampleRate/(double)fftlen);
   outtextxy(FRX,FRY,ach);
   sprintf(ach,"%s window function",window_name[windfunc]);
   outtextxy(WFX,WFY,ach);
   sprintf(ach,"Gain of %d dB per octave",gain);
   outtextxy(DGX,DGY,ach);
   sprintf(ach,"(0 dB gain at %ld Hz)",ref_freq);
   outtextxy(RFX,RFY,ach);
}