rlm@23: #include <stdio.h>
rlm@23: #include <stdlib.h>
rlm@23: #include <math.h>
rlm@23: #include <assert.h>
rlm@23: #include <string.h>
rlm@23: #include <math.h>
rlm@23: #include "SndfileWavUtil.h"
rlm@23: 
rlm@23: const  int Taps = 9;
rlm@23: const  int Points = 8;
rlm@23: const  int log_Points = 3;
rlm@23: static int newest;
rlm@23: static short int x[Taps] = {0,0,0,0,0,0,0,0,0};
rlm@23: static short int fir(short int);
rlm@23: 
rlm@23: // foreward declaration. This is defined in DFT.cpp
rlm@23: int DFT(int dir,int m,double *x1,double *y1);
rlm@23: 
rlm@23: 
rlm@23: static int x_to_y(int x, int y)
rlm@23: {
rlm@23:   assert(y>0);
rlm@23:   int rv = x;
rlm@23:   while(y-->1)
rlm@23:     x *= x;
rlm@23:   return rv;
rlm@23: }
rlm@23: 
rlm@23: // function implemented in LinearWindowingFunction.bsv
rlm@23: static double windowing_fn(int index)
rlm@23: {
rlm@23:   double pi = 3.14159265;
rlm@23:   double i = (double)index;
rlm@23:   return sin(0.5*pi* (sin((i + 0.5)/Points * pi) * sin((i + 0.5)/Points * pi)));
rlm@23: }
rlm@23: 
rlm@23: // lifted from the bluespec code
rlm@23: static double h[Taps] =  {-0.0124,
rlm@23: 			  0.0,
rlm@23: 			  -0.0133,
rlm@23: 			  0.0,
rlm@23: 			  0.8181,
rlm@23: 			  0.0,
rlm@23: 			  -0.0133,
rlm@23: 			  0.0,
rlm@23: 			  -0.0124};
rlm@23: 
rlm@23: 
rlm@23: short int
rlm@23: fir(short int sample)
rlm@23: {
rlm@23:   x[newest] = sample; 
rlm@23:   double y = 0;
rlm@23:   for (int k = 0; k < Taps; k++) { 
rlm@23:     short int bits = x[(newest+k)%Taps];
rlm@23:     double x = (double)bits;
rlm@23:     y += h[k] * x; 
rlm@23:   }
rlm@23:   newest = (newest == 0) ? (Taps-1) : (newest-1);
rlm@23:   short int rv = ((short int)y)&0x0000FFFF;
rlm@23:   return rv;
rlm@23: }
rlm@23: 
rlm@23: 
rlm@23: int
rlm@23: main (int argc, char * argv [])
rlm@23: {
rlm@23:   const char* inputWavFileName;
rlm@23:   const char* outputWavFileName;
rlm@23: 
rlm@23:   FILE *inputPcmFile;
rlm@23:   FILE *firPcmFile;
rlm@23:   FILE *ifftPcmFile;
rlm@23:   FILE *outputPcmFile;
rlm@23: 
rlm@23:   short int sample;
rlm@23:   short int samples[Points];
rlm@23:   const unsigned int short_sz = sizeof(short int);
rlm@23: 
rlm@23:   inputWavFileName = argv[1];
rlm@23:   outputWavFileName = argv[2];
rlm@23: 
rlm@23:   // Convert input wav to pcm
rlm@23:   generate_pcm(inputWavFileName, "input.pcm");
rlm@23: 
rlm@23:   inputPcmFile  = fopen("input.pcm", "r");
rlm@23:   firPcmFile    = fopen("fir.pcm", "w");
rlm@23: 
rlm@23:   assert(inputPcmFile);
rlm@23:   assert(firPcmFile);
rlm@23: 
rlm@23:   newest = 0;
rlm@23:   memset(x, 0, sizeof(x));
rlm@23: 
rlm@23:   while(fread(&sample, short_sz, 1, inputPcmFile)) {
rlm@23:     sample = fir(sample);
rlm@23:     assert(fwrite(&sample,short_sz, 1, firPcmFile));
rlm@23:   }
rlm@23: 
rlm@23:   fclose(inputPcmFile);
rlm@23:   fclose(firPcmFile);
rlm@23: 
rlm@23:   firPcmFile  = fopen("fir.pcm", "r");
rlm@23:   ifftPcmFile = fopen("ifft.pcm", "w");
rlm@23: 
rlm@23:   assert(firPcmFile);
rlm@23:   assert(ifftPcmFile);
rlm@23: 
rlm@23:   int read = 0;
rlm@23: 
rlm@23:   // read in half a frame
rlm@23:   for(int i = 0; i < Points/2; i++)
rlm@23:     read += fread(&samples[(Points/2)+i], short_sz, 1, firPcmFile);
rlm@23: 
rlm@23:   if(read > 0)
rlm@23:     assert(read==Points/2);
rlm@23: 
rlm@23:   // we will do an 'Points' point fft, and then undo it
rlm@23:   while(true) {
rlm@23: 
rlm@23:     read = 0;
rlm@23: 
rlm@23:     // shift last half frame to head
rlm@23:     for(int i = 0; i < Points/2; i++)
rlm@23:       samples[i] = samples[(Points/2)+i];
rlm@23: 
rlm@23:     // read in another half frame
rlm@23:     for(int i = 0; i < Points/2; i++)
rlm@23:       read += fread(&samples[(Points/2)+i], short_sz, 1, firPcmFile);
rlm@23: 
rlm@23:     if(read == 0)
rlm@23:       break;
rlm@23: 
rlm@23:     // pad out the 
rlm@23:     if(read < Points/2){
rlm@23:       for(int i = read; i < Points/2; i++)
rlm@23: 	samples[Points/2+i] = 0;
rlm@23:     }
rlm@23: 
rlm@23:     double dsamples[Points];
rlm@23:     double dimag[Points];
rlm@23: 
rlm@23:     // this shift is performed in the Bluespec
rlm@23:     for(int i = 0; i < Points; i++){
rlm@23:       dsamples[i] = (double)(samples[i]>>log_Points);
rlm@23:       dimag[i] = 0.0;
rlm@23:     }
rlm@23: 
rlm@23: 
rlm@23:     //for(int i = 0; i < Points; i++)
rlm@23:     //  printf("%d ", (int)dsamples[i]);
rlm@23:     //printf("\n");
rlm@23: 
rlm@23:     DFT(1,Points,dsamples,dimag);
rlm@23: 
rlm@23:     // shift up by 'shift' (in)
rlm@23:     int shift = 2;
rlm@23: 
rlm@23:     for(int i = 0; i < Points/2; i++){
rlm@23: 
rlm@23:       dsamples[Points/2-1-i] = dsamples[Points/2-1-i-shift];
rlm@23:       dimag[Points/2-1-i] = dimag[Points/2-1-i-shift];
rlm@23: 
rlm@23: 
rlm@23:       dsamples[Points/2+i] = dsamples[Points/2+i+shift];
rlm@23:       dimag[Points/2+i] = dimag[Points/2+i+shift];
rlm@23: 
rlm@23:     }
rlm@23: 
rlm@23:     // fill in the lower points
rlm@23:     for(int i = 0; i < shift; i++){
rlm@23: 
rlm@23:       dsamples[i] = 0.0;
rlm@23:       dimag[i] = 0.0;
rlm@23: 
rlm@23:       dsamples[Points-1-i] = 0.0;
rlm@23:       dimag[Points-1-i] = 0.0;
rlm@23:     }
rlm@23: 
rlm@23:     DFT(-1,Points,dsamples,dimag);
rlm@23: 
rlm@23:     short int ifftResult[Points];
rlm@23:     for(int i = 0; i < Points; i++)
rlm@23:       ifftResult[i] = (int)dsamples[i];
rlm@23: 
rlm@23:     // strip off the padding
rlm@23:     int write = (read < Points/2) ? read+(Points/2) : Points;
rlm@23:     for(int i = 0; i < write; i++)
rlm@23:       assert(fwrite(&ifftResult[i],short_sz,1,ifftPcmFile));
rlm@23: 
rlm@23:     // break out if we're at the end
rlm@23:     if(read < Points/2){
rlm@23:        break;
rlm@23:     }
rlm@23:   }
rlm@23: 
rlm@23:   fclose(firPcmFile);
rlm@23:   fclose(ifftPcmFile);
rlm@23: 
rlm@23:   // do the windowing
rlm@23: 
rlm@23:   ifftPcmFile   = fopen("ifft.pcm", "r");
rlm@23:   outputPcmFile = fopen("output.pcm", "w");
rlm@23: 
rlm@23:   assert(ifftPcmFile);
rlm@23:   assert(outputPcmFile);
rlm@23: 
rlm@23:   short int samplesA[Points];
rlm@23:   short int samplesB[Points];
rlm@23: 
rlm@23:   bool valid_stream = true;
rlm@23:   // write out the first half frame
rlm@23:   for(int i = 0; i < Points/2; i++){
rlm@23:     valid_stream &= fread(samplesA, short_sz, 1,ifftPcmFile);
rlm@23:     valid_stream &= fwrite(samplesA, short_sz,1,outputPcmFile);
rlm@23:   }
rlm@23: 
rlm@23:   while(valid_stream){
rlm@23: 
rlm@23:     for(int i = 0; i < Points/2; i++)
rlm@23:       valid_stream &= fread(&samplesA[i], short_sz, 1,ifftPcmFile);
rlm@23: 
rlm@23:     for(int i = 0; i < Points; i++)
rlm@23:       valid_stream &= fread(&samplesB[i], short_sz, 1,ifftPcmFile);
rlm@23: 
rlm@23:     for(int i = 0; i < Points/2; i++)
rlm@23:       valid_stream &= fread(&samplesA[(Points/2)+i], short_sz, 1,ifftPcmFile); 
rlm@23: 
rlm@23: 
rlm@23:     // this isn't quite right
rlm@23:     if(!valid_stream)
rlm@23:       break;
rlm@23: 
rlm@23:     for(int i = 0; i < Points; i++){
rlm@23:       double window = windowing_fn(i);
rlm@23:       short int rv = (short int)(((double)samplesA[i])*(1.0-window) + window*((double)samplesB[i]));
rlm@23:       assert(fwrite(&rv,short_sz,1,outputPcmFile));    
rlm@23:     }
rlm@23: 
rlm@23:   }
rlm@23: 
rlm@23:   fclose(ifftPcmFile);
rlm@23:   fclose(outputPcmFile);
rlm@23: 
rlm@23:   generate_wav("output.pcm", inputWavFileName, outputWavFileName);
rlm@23: 
rlm@23: }
rlm@23: