#include "nrutil.h"
#include<stdio.h>
#include <malloc.h>
#include <math.h> 

#define HCONST (60.0) /* Constant used in plotting the spectrum */


/* ============ Generate LPC and/or FFT plots ================= */

void MAKE_PLOTS(a,s,gain,fftSpec,fp,fpgnu,fspec,N,NP,F_LEN)
float *a,*s,gain;
int fftSpec,N,NP,F_LEN;
FILE *fp,*fpgnu,*fspec;
{
int i;
static int from=0;
static int cc=1;
float *freq;
void FREQ_RESPONSE();
void FFT_SPECTRUM();


   freq = vector (1,N/2+1);

    FREQ_RESPONSE(a,gain,freq,N,NP);

    for (i=1;i<=N/2;++i)
      fprintf(fspec,"%f\n",freq[i]); 

    if (fftSpec){ /*-- print the FFT spectrum as well --- */
      FFT_SPECTRUM(s,fp,N,F_LEN);
      fprintf(fpgnu,"set title 'Frame  %2d'\n",cc);
      fprintf(fpgnu,"plot [%d:%d][] 'spec.dat' w l,'fft.dat' w l\n",from,
	      from+N/2-1); }
    else
      fprintf(fpgnu,"plot [%d:%d][] 'spec.dat' w l\n",from,from+N/2-1);
   

    fprintf(fpgnu,"pause -1\n"); 
    from += N/2;


 free_vector(freq,1,N/2+1);
cc++;
}


/* ================================================== */
void FREQ_RESPONSE(a,gain,freq,N,NP)
float *a,*freq,gain;
int N,NP;
{
float *sx,bin;
int i,k,first,start,stop,now;
float real,imag,dB,magnitude,inverse,temp;
void FFT();


sx = vector(1,2*N+2); 

/* ----------- read in the lpc coefficients ----------------- */
/**/
k=3;
sx[1]=1.0;
sx[2]=0.0;
for(i=2;i<=N; ++i)
  {
   if(i<=NP+1){
     sx[k]=a[i-1];
     sx[k+1]=0.0; 
   }
   else{
     sx[k]=0.0;
     sx[k+1]=0.0;
   }
  k +=2;
}    


 FFT(sx,N,1);

/* -------- Compute frequency response of 1/H(z) -------------- */
 k=1;
 for(i=1;i<=N/2;++i){
   real=  sx[k];
   imag = sx[k+1];
   magnitude=sqrt(real*real+imag*imag); 
   inverse=gain/magnitude;
   /* dB = 10.0*log10(inverse*gain);  */
   freq[i]= inverse;
   k+=2;
 }

 free_vector(sx,1,2*N+2);

}
/************************************************************* */
void FFT_SPECTRUM(sp,fp1,N,F_LEN)
float *sp;
FILE *fp1;
int N,F_LEN;
{
FILE *fp2;
float *sin,*freq,bin;
int i,k,first,start,stop,now;
float real,imag,dB,magnitude,inverse,temp;
void FFT();


 sin  = vector(1,2*N+2);
 freq = vector(1,N/2+1);

/* ----------- read in speech -------------- */

k=1;
for(i=1;i<=N; ++i)
  {
    if (i<=F_LEN) {
     sin[k]  =sp[i-1];
     sin[k+1]=0.0;
   }
    else {
     sin[k]=0.0;
     sin[k+1]=0.0;
   }

  k +=2;
}    
 

 FFT(sin,N,1);

/* -------- Compute frequency response -------------- */

 k=1;
 for(i=1;i<=N/2;++i){
   real=sin[k];
   imag = sin[k+1];
   magnitude=real*real+imag*imag;
   dB = 10.0*log10(magnitude/HCONST); 
   fprintf(fp1,"%f\n",dB);
   k+=2;
 }

}