Mercurial > pygar

     1 #include <stdio.h>

     2 #include <stdlib.h>

     3 #include <math.h>

     4 #include <assert.h>

     5 #include <string.h>

     6 #include "SndfileWavUtil.h"

     7 

     8 const  int Taps = 9;

     9 const  int Points = 8;

    10 const  int log_Points = 3;

    11 static int newest;

    12 static short int x[Taps] = {0,0,0,0,0,0,0,0,0};

    13 static short int fir(short int);

    14 

    15 // foreward declaration. This is defined in DFT.cpp

    16 int DFT(int dir,int m,double *x1,double *y1);

    17 

    18 // lifted from the bluespec code

    19 static double h[Taps] =  {-0.0124,

    20 			  0.0,

    21 			  -0.0133,

    22 			  0.0,

    23 			  0.8181,

    24 			  0.0,

    25 			  -0.0133,

    26 			  0.0,

    27 			  -0.0124};

    28 

    29 int

    30 main (int argc, char * argv [])

    31 {

    32   const char* inputWavFileName;

    33   const char* outputWavFileName;

    34 

    35   FILE *inputPcmFile;

    36   FILE *firPcmFile;

    37   FILE *outputPcmFile;

    38 

    39   short int sample;

    40   short int samples[Points];

    41 

    42   inputWavFileName = argv[1];

    43   outputWavFileName = argv[2];

    44 

    45   // Convert input wav to pcm

    46   generate_pcm(inputWavFileName, "input.pcm");

    47 

    48   inputPcmFile  = fopen("input.pcm", "r");

    49   firPcmFile    = fopen("fir.pcm", "w");

    50 

    51   assert(inputPcmFile);

    52   assert(firPcmFile);

    53 

    54   newest = 0;

    55   memset(x, 0, sizeof(x));

    56 

    57   while(fread(&sample, 2, 1, inputPcmFile)) {

    58     sample = fir(sample);

    59     assert(fwrite(&sample,2,1,firPcmFile));

    60   }

    61 

    62   fclose(inputPcmFile);

    63   fclose(firPcmFile);

    64 

    65   firPcmFile    = fopen("fir.pcm", "r");

    66   outputPcmFile = fopen("output.pcm", "w");

    67 

    68   assert(firPcmFile);

    69   assert(outputPcmFile);

    70 

    71   // we will do an Points point fft, and then undo it

    72   while(true) {

    73 

    74     int read = fread(samples, 2,Points,firPcmFile);

    75 

    76     if(read == 0)

    77       break;

    78 

    79     // pad out the 

    80     if(read<Points){

    81       for(int i = read; i < Points; i++)

    82 	samples[i] = 0;

    83     }

    84 

    85     double dsamples[Points];

    86     double dimag[Points];

    87 

    88     // This shift is performed in the HW

    89     for(int i = 0; i < Points; i++){

    90       dsamples[i] = (double)(samples[i]>>log_Points);

    91       dimag[i] = 0.0;

    92     }

    93 

    94     DFT(1,Points,dsamples,dimag);

    95     DFT(-1,Points,dsamples,dimag);

    96 

    97     for(int i = 0; i < Points; i++){

    98       samples[i] = (int)dsamples[i];

    99     }

   100 

   101     if(read < Points){

   102       assert(fwrite(samples,2,read,outputPcmFile)==read);

   103       break;

   104     } else {

   105       assert(fwrite(samples,2,Points,outputPcmFile)==Points);

   106     }

   107   }

   108 

   109   fclose(firPcmFile);

   110   fclose(outputPcmFile);

   111 

   112   generate_wav("output.pcm", inputWavFileName, outputWavFileName);

   113 

   114 }

   115 

   116 short int

   117 fir(short int sample)

   118 {

   119 

   120   x[newest] = sample; 

   121   double y = 0;

   122 

   123   for (int k = 0; k < Taps; k++) { 

   124 

   125     short int bits = x[(newest+k)%Taps];

   126     double x = (double)bits;

   127     y += h[k] * x; 

   128 

   129   }

   130 

   131   newest = (newest == 0) ? (Taps-1) : (newest-1);

   132 

   133   short int rv = ((short int)y)&0x0000FFFF;

   134   return rv;

   135 }

   136
author	rlm
date	Tue, 11 May 2010 23:23:21 -0400
parents	90197e3375e2
children