/* GENKERN.C  4/91 brown@nadia
        usage: genkern -flags > outfile of floats
  genkern > del ;  ; fft -fR 256 256 < del > kernR.256.60
  "kernR.256.60" is read by cqtrans to take the const Q trans.

/* generate kernels which are transforms of    w(t)exp(-j wo t)   where w(t)
is has width 20ms for freqs up to to 1667 Hz  and then corresp to qtrron
 resolution; SR = 11025 ;
Freq wo from F3 to Nyquist with freqs corresp to  = tempered scale.
NB: tuning of the channels dep on MINHZ or srate and can be quite critical.
*/
#include <stdio.h>
#include <math.h>
#define Nwind 256
#define MINHZ 174.6  /* F3 = 7.05  midi=53 */
#define TWOPI  6.2831854 
#define PI  3.1415927 
#define NMAXFREQS 150


main(argc,argv)
 int argc;
 char **argv;
{  /* begin main */
        char *cp;
        int flag, debug=0 , hamming = 1, dbout=0;
	int k, comp, n, topval=0 , varwind=1, lowval=0, nozeros=0,wndsz;
	int   setzero=0, CnstWind = Nwind;
	int windsiz[NMAXFREQS];
	float srate = 11025. ;
	float ETF[NMAXFREQS];
	float kernel[20000];
	float windshape[Nwind];
        float theta, two1twlfth = pow(2.,(1./12.)), *kernelp = kernel;
	float two1tfrth = pow(2.,(1./24.));
	float mindigfreq, alpha, factor, twopidws = TWOPI/Nwind;
	float osix = two1twlfth - 1. , oneMalpha, phi;
	float othree = two1tfrth - 1., res;

        argc--; argv++;                         /* align onto command args */
        while ((cp = *argv) && *cp++ =='-' && (flag = *cp++))
        {
                switch (flag){
                        case 'd': sscanf(cp, "%f", &debug);
                                break;
                        case 'h': sscanf(cp, "%f", &hamming);
                                break;
                        case 't': sscanf(cp, "%d", &topval);
                                break;               /* highest kernel freq */
                        case 'l': sscanf(cp, "%d", &lowval);
                                break;                /* lowest kernel freq */
                        default: fprintf(stderr, "unknown option\n");
                                break;
                }
                argc--; argv++;
        }
	if(hamming) 
		 alpha=(double)25./(double)46.;
	oneMalpha = 1. - alpha ;
	res =  othree ;
	fprintf(stderr,"genkern: %s window, %s out.\n",
		(hamming) ? "hamming":"rect", (dbout) ? "db":"magnitude");

/******  generate = tempered freqs and windsiz[k] ******************/
	mindigfreq = TWOPI * MINHZ /srate ; /*lowest dig freq -> F3 */
	for (k = 0; k < NMAXFREQS ; ++k){
	   ETF[k] = mindigfreq * pow(two1twlfth, (float)k); 
	                     /* windsiz[k]=srate/.06f = TWOPI/(ETF[k]*res */
	   windsiz[k] = (int)(TWOPI/(ETF[k] * res)) ;
 	   if (debug)
	     fprintf(stderr,"k = %d  digfreq = %.3f freq=%.0f windsiz=%d\n", k, ETF[k], srate* ETF[k]/TWOPI, windsiz[k] ); 

	   if(ETF[k] > PI )
	     break;
	}  /********  remember kth value is over the Nyquist *******/

/************  calc topval minus lowval complex kernels *******/
	if (!topval)
	  topval = k; /* option so few enough frames for graph viewable */
	fprintf(stderr,"Lowest comp %d; highest comp %d\n",lowval, topval);
	for (comp = lowval; comp < topval; ++comp){
	   for (n = 0; n < Nwind; ++n){ 
	      setzero = 0;
	      if (hamming){
		  wndsz=(windsiz[comp]<CnstWind ? windsiz[comp]-1 :CnstWind-1);
		  nozeros = CnstWind - wndsz;
		  if((n < nozeros/2) || (n > ((CnstWind-1)- nozeros/2) )){
		     phi = 0.;
		     setzero = 1;
		                /* to get zeros in the window outside Hamming*/
		  }
		  else
		    phi = (n - nozeros/2)*TWOPI/wndsz ;
		                       /* displace n so hamming is centered */
	      }
	      if(setzero)
		windshape[n] = 0.;  /* windshape doesn't need to be an array*/
	      else
		windshape[n] =(!hamming)? 1. : alpha- (oneMalpha * cos(phi));
	      theta = ETF[comp] * (n- CnstWind/2) ;        /* center the fn */
	      *kernelp  =(windshape[n] * cos(theta)) /wndsz; /* wndsz fr norm*/
	      ++kernelp;
	   }/* end n loop */
	}
	write(1, kernel, (topval - lowval)*Nwind*4); /* just one frm pr com*/

} /* end main */

merde(s)
 char *s;
{
        fprintf(stderr,"MERDE, ALORS!   %s\n",s);
        exit(1);
}