1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/tools/audio_processor_test/transform/checker/checker.cpp	Wed Apr 28 08:19:09 2010 -0400
     1.3 @@ -0,0 +1,219 @@
     1.4 +#include <stdio.h>
     1.5 +#include <stdlib.h>
     1.6 +#include <math.h>
     1.7 +#include <assert.h>
     1.8 +#include <string.h>
     1.9 +#include "SndfileWavUtil.h"
    1.10 +
    1.11 +const  int Taps = 9;
    1.12 +const  int Points = 8;
    1.13 +const  int log_Points = 3;
    1.14 +static int newest;
    1.15 +static short int x[Taps] = {0,0,0,0,0,0,0,0,0};
    1.16 +static short int fir(short int);
    1.17 +
    1.18 +// foreward declaration. This is defined in DFT.cpp
    1.19 +int DFT(int dir,int m,double *x1,double *y1);
    1.20 +
    1.21 +
    1.22 +static int x_to_y(int x, int y)
    1.23 +{
    1.24 +  assert(y>0);
    1.25 +  int rv = x;
    1.26 +  while(y-->1)
    1.27 +    x *= x;
    1.28 +  return rv;
    1.29 +}
    1.30 +
    1.31 +// function implemented in LinearWindowingFunction.bsv
    1.32 +static double windowing_fn(int index)
    1.33 +{
    1.34 +  double i = (double)index;
    1.35 +  double divisor = Points/2; //(double) (x_to_y(2,log_Points-1));
    1.36 +  double rv = (index < Points/2) ? (i / divisor) : ((((double)Points)-i) / divisor);
    1.37 +  return rv;
    1.38 +}
    1.39 +
    1.40 +// lifted from the bluespec code
    1.41 +static double h[Taps] =  {-0.0124,
    1.42 +			  0.0,
    1.43 +			  -0.0133,
    1.44 +			  0.0,
    1.45 +			  0.8181,
    1.46 +			  0.0,
    1.47 +			  -0.0133,
    1.48 +			  0.0,
    1.49 +			  -0.0124};
    1.50 +
    1.51 +
    1.52 +short int
    1.53 +fir(short int sample)
    1.54 +{
    1.55 +  x[newest] = sample; 
    1.56 +  double y = 0;
    1.57 +  for (int k = 0; k < Taps; k++) { 
    1.58 +    short int bits = x[(newest+k)%Taps];
    1.59 +    double x = (double)bits;
    1.60 +    y += h[k] * x; 
    1.61 +  }
    1.62 +  newest = (newest == 0) ? (Taps-1) : (newest-1);
    1.63 +  short int rv = ((short int)y)&0x0000FFFF;
    1.64 +  return rv;
    1.65 +}
    1.66 +
    1.67 +
    1.68 +int
    1.69 +main (int argc, char * argv [])
    1.70 +{
    1.71 +  const char* inputWavFileName;
    1.72 +  const char* outputWavFileName;
    1.73 +
    1.74 +  FILE *inputPcmFile;
    1.75 +  FILE *firPcmFile;
    1.76 +  FILE *ifftPcmFile;
    1.77 +  FILE *outputPcmFile;
    1.78 +
    1.79 +  short int sample;
    1.80 +  short int samples[Points];
    1.81 +  const unsigned int short_sz = sizeof(short int);
    1.82 +
    1.83 +  inputWavFileName = argv[1];
    1.84 +  outputWavFileName = argv[2];
    1.85 +
    1.86 +  // Convert input wav to pcm
    1.87 +  generate_pcm(inputWavFileName, "input.pcm");
    1.88 +
    1.89 +  inputPcmFile  = fopen("input.pcm", "r");
    1.90 +  firPcmFile    = fopen("fir.pcm", "w");
    1.91 +
    1.92 +  assert(inputPcmFile);
    1.93 +  assert(firPcmFile);
    1.94 +
    1.95 +  newest = 0;
    1.96 +  memset(x, 0, sizeof(x));
    1.97 +
    1.98 +  while(fread(&sample, short_sz, 1, inputPcmFile)) {
    1.99 +    sample = fir(sample);
   1.100 +    assert(fwrite(&sample,short_sz, 1, firPcmFile));
   1.101 +  }
   1.102 +
   1.103 +  fclose(inputPcmFile);
   1.104 +  fclose(firPcmFile);
   1.105 +
   1.106 +  firPcmFile  = fopen("fir.pcm", "r");
   1.107 +  ifftPcmFile = fopen("ifft.pcm", "w");
   1.108 +
   1.109 +  assert(firPcmFile);
   1.110 +  assert(ifftPcmFile);
   1.111 +
   1.112 +  int read = 0;
   1.113 +
   1.114 +  // read in half a frame
   1.115 +  for(int i = 0; i < Points/2; i++)
   1.116 +    read += fread(&samples[(Points/2)+i], short_sz, 1, firPcmFile);
   1.117 +
   1.118 +  if(read > 0)
   1.119 +    assert(read==Points/2);
   1.120 +
   1.121 +  // we will do an 'Points' point fft, and then undo it
   1.122 +  while(true) {
   1.123 +
   1.124 +    read = 0;
   1.125 +
   1.126 +    // shift last half frame to head
   1.127 +    for(int i = 0; i < Points/2; i++)
   1.128 +      samples[i] = samples[(Points/2)+i];
   1.129 +
   1.130 +    // read in another half frame
   1.131 +    for(int i = 0; i < Points/2; i++)
   1.132 +      read += fread(&samples[(Points/2)+i], short_sz, 1, firPcmFile);
   1.133 +
   1.134 +    if(read == 0)
   1.135 +      break;
   1.136 +
   1.137 +    // pad out the 
   1.138 +    if(read < Points/2){
   1.139 +      for(int i = read; i < Points/2; i++)
   1.140 +	samples[Points/2+i] = 0;
   1.141 +    }
   1.142 +
   1.143 +    double dsamples[Points];
   1.144 +    double dimag[Points];
   1.145 +
   1.146 +    // this shift is performed in the Bluespec
   1.147 +    for(int i = 0; i < Points; i++){
   1.148 +      dsamples[i] = (double)(samples[i]>>log_Points);
   1.149 +      dimag[i] = 0.0;
   1.150 +    }
   1.151 +
   1.152 +
   1.153 +    DFT(1,Points,dsamples,dimag);
   1.154 +    DFT(-1,Points,dsamples,dimag);
   1.155 +
   1.156 +    short int ifftResult[Points];
   1.157 +    for(int i = 0; i < Points; i++)
   1.158 +      ifftResult[i] = (int)dsamples[i]; //)>>log_Points;
   1.159 +
   1.160 +    // strip off the padding
   1.161 +    int write = (read < Points/2) ? read+(Points/2) : Points;
   1.162 +    for(int i = 0; i < write; i++)
   1.163 +      assert(fwrite(&ifftResult[i],short_sz,1,ifftPcmFile));
   1.164 +
   1.165 +    // break out if we're at the end
   1.166 +    if(read < Points/2){
   1.167 +       break;
   1.168 +    }
   1.169 +  }
   1.170 +
   1.171 +  fclose(firPcmFile);
   1.172 +  fclose(ifftPcmFile);
   1.173 +
   1.174 +  // do the windowing
   1.175 +
   1.176 +  ifftPcmFile   = fopen("ifft.pcm", "r");
   1.177 +  outputPcmFile = fopen("output.pcm", "w");
   1.178 +
   1.179 +  assert(ifftPcmFile);
   1.180 +  assert(outputPcmFile);
   1.181 +
   1.182 +  short int samplesA[Points];
   1.183 +  short int samplesB[Points];
   1.184 +
   1.185 +  bool valid_stream = true;
   1.186 +  // write out the first half frame
   1.187 +  for(int i = 0; i < Points/2; i++){
   1.188 +    valid_stream &= fread(samplesA, short_sz, 1,ifftPcmFile);
   1.189 +    valid_stream &= fwrite(samplesA, short_sz,1,outputPcmFile);
   1.190 +  }
   1.191 +
   1.192 +  while(valid_stream){
   1.193 +
   1.194 +    for(int i = 0; i < Points/2; i++)
   1.195 +      valid_stream &= fread(&samplesA[i], short_sz, 1,ifftPcmFile);
   1.196 +
   1.197 +    for(int i = 0; i < Points; i++)
   1.198 +      valid_stream &= fread(&samplesB[i], short_sz, 1,ifftPcmFile);
   1.199 +
   1.200 +    for(int i = 0; i < Points/2; i++)
   1.201 +      valid_stream &= fread(&samplesA[(Points/2)+i], short_sz, 1,ifftPcmFile); 
   1.202 +
   1.203 +
   1.204 +    // this isn't quite right
   1.205 +    if(!valid_stream)
   1.206 +      break;
   1.207 +
   1.208 +    for(int i = 0; i < Points; i++){
   1.209 +      double window = windowing_fn(i);
   1.210 +      short int rv = samplesA[i]*(1-window) + window*samplesB[i];
   1.211 +      assert(fwrite(&rv,short_sz,1,outputPcmFile));    
   1.212 +    }
   1.213 +
   1.214 +  }
   1.215 +
   1.216 +  fclose(ifftPcmFile);
   1.217 +  fclose(outputPcmFile);
   1.218 +
   1.219 +  generate_wav("output.pcm", inputWavFileName, outputWavFileName);
   1.220 +
   1.221 +}
   1.222 +