Mercurial > pygar

annotate tools/audio_processor_test/cosine_window/checker/checker.cpp @ 59:92041177735c pygar svn.60

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[svn r60] LUTRAMing
author punk
date Mon, 10 May 2010 13:50:40 -0400
parents 90197e3375e2
children
rev   line source
rlm@23 1 #include <stdio.h>
rlm@23 2 #include <stdlib.h>
rlm@23 3 #include <math.h>
rlm@23 4 #include <assert.h>
rlm@23 5 #include <string.h>
rlm@23 6 #include <math.h>
rlm@23 7 #include "SndfileWavUtil.h"
rlm@23 8
rlm@23 9 const int Taps = 9;
rlm@23 10 const int Points = 8;
rlm@23 11 const int log_Points = 3;
rlm@23 12 static int newest;
rlm@23 13 static short int x[Taps] = {0,0,0,0,0,0,0,0,0};
rlm@23 14 static short int fir(short int);
rlm@23 15
rlm@23 16 // foreward declaration. This is defined in DFT.cpp
rlm@23 17 int DFT(int dir,int m,double *x1,double *y1);
rlm@23 18
rlm@23 19
rlm@23 20 static int x_to_y(int x, int y)
rlm@23 21 {
rlm@23 22 assert(y>0);
rlm@23 23 int rv = x;
rlm@23 24 while(y-->1)
rlm@23 25 x *= x;
rlm@23 26 return rv;
rlm@23 27 }
rlm@23 28
rlm@23 29 // function implemented in LinearWindowingFunction.bsv
rlm@23 30 static double windowing_fn(int index)
rlm@23 31 {
rlm@23 32 double pi = 3.14159265;
rlm@23 33 double i = (double)index;
rlm@23 34 return sin(0.5*pi* (sin((i + 0.5)/Points * pi) * sin((i + 0.5)/Points * pi)));
rlm@23 35 }
rlm@23 36
rlm@23 37 // lifted from the bluespec code
rlm@23 38 static double h[Taps] = {-0.0124,
rlm@23 39 0.0,
rlm@23 40 -0.0133,
rlm@23 41 0.0,
rlm@23 42 0.8181,
rlm@23 43 0.0,
rlm@23 44 -0.0133,
rlm@23 45 0.0,
rlm@23 46 -0.0124};
rlm@23 47
rlm@23 48
rlm@23 49 short int
rlm@23 50 fir(short int sample)
rlm@23 51 {
rlm@23 52 x[newest] = sample;
rlm@23 53 double y = 0;
rlm@23 54 for (int k = 0; k < Taps; k++) {
rlm@23 55 short int bits = x[(newest+k)%Taps];
rlm@23 56 double x = (double)bits;
rlm@23 57 y += h[k] * x;
rlm@23 58 }
rlm@23 59 newest = (newest == 0) ? (Taps-1) : (newest-1);
rlm@23 60 short int rv = ((short int)y)&0x0000FFFF;
rlm@23 61 return rv;
rlm@23 62 }
rlm@23 63
rlm@23 64
rlm@23 65 int
rlm@23 66 main (int argc, char * argv [])
rlm@23 67 {
rlm@23 68 const char* inputWavFileName;
rlm@23 69 const char* outputWavFileName;
rlm@23 70
rlm@23 71 FILE *inputPcmFile;
rlm@23 72 FILE *firPcmFile;
rlm@23 73 FILE *ifftPcmFile;
rlm@23 74 FILE *outputPcmFile;
rlm@23 75
rlm@23 76 short int sample;
rlm@23 77 short int samples[Points];
rlm@23 78 const unsigned int short_sz = sizeof(short int);
rlm@23 79
rlm@23 80 inputWavFileName = argv[1];
rlm@23 81 outputWavFileName = argv[2];
rlm@23 82
rlm@23 83 // Convert input wav to pcm
rlm@23 84 generate_pcm(inputWavFileName, "input.pcm");
rlm@23 85
rlm@23 86 inputPcmFile = fopen("input.pcm", "r");
rlm@23 87 firPcmFile = fopen("fir.pcm", "w");
rlm@23 88
rlm@23 89 assert(inputPcmFile);
rlm@23 90 assert(firPcmFile);
rlm@23 91
rlm@23 92 newest = 0;
rlm@23 93 memset(x, 0, sizeof(x));
rlm@23 94
rlm@23 95 while(fread(&sample, short_sz, 1, inputPcmFile)) {
rlm@23 96 sample = fir(sample);
rlm@23 97 assert(fwrite(&sample,short_sz, 1, firPcmFile));
rlm@23 98 }
rlm@23 99
rlm@23 100 fclose(inputPcmFile);
rlm@23 101 fclose(firPcmFile);
rlm@23 102
rlm@23 103 firPcmFile = fopen("fir.pcm", "r");
rlm@23 104 ifftPcmFile = fopen("ifft.pcm", "w");
rlm@23 105
rlm@23 106 assert(firPcmFile);
rlm@23 107 assert(ifftPcmFile);
rlm@23 108
rlm@23 109 int read = 0;
rlm@23 110
rlm@23 111 // read in half a frame
rlm@23 112 for(int i = 0; i < Points/2; i++)
rlm@23 113 read += fread(&samples[(Points/2)+i], short_sz, 1, firPcmFile);
rlm@23 114
rlm@23 115 if(read > 0)
rlm@23 116 assert(read==Points/2);
rlm@23 117
rlm@23 118 // we will do an 'Points' point fft, and then undo it
rlm@23 119 while(true) {
rlm@23 120
rlm@23 121 read = 0;
rlm@23 122
rlm@23 123 // shift last half frame to head
rlm@23 124 for(int i = 0; i < Points/2; i++)
rlm@23 125 samples[i] = samples[(Points/2)+i];
rlm@23 126
rlm@23 127 // read in another half frame
rlm@23 128 for(int i = 0; i < Points/2; i++)
rlm@23 129 read += fread(&samples[(Points/2)+i], short_sz, 1, firPcmFile);
rlm@23 130
rlm@23 131 if(read == 0)
rlm@23 132 break;
rlm@23 133
rlm@23 134 // pad out the
rlm@23 135 if(read < Points/2){
rlm@23 136 for(int i = read; i < Points/2; i++)
rlm@23 137 samples[Points/2+i] = 0;
rlm@23 138 }
rlm@23 139
rlm@23 140 double dsamples[Points];
rlm@23 141 double dimag[Points];
rlm@23 142
rlm@23 143 // this shift is performed in the Bluespec
rlm@23 144 for(int i = 0; i < Points; i++){
rlm@23 145 dsamples[i] = (double)(samples[i]>>log_Points);
rlm@23 146 dimag[i] = 0.0;
rlm@23 147 }
rlm@23 148
rlm@23 149
rlm@23 150 //for(int i = 0; i < Points; i++)
rlm@23 151 // printf("%d ", (int)dsamples[i]);
rlm@23 152 //printf("\n");
rlm@23 153
rlm@23 154 DFT(1,Points,dsamples,dimag);
rlm@23 155
rlm@23 156 // shift up by 'shift' (in)
rlm@23 157 int shift = 2;
rlm@23 158
rlm@23 159 for(int i = 0; i < Points/2; i++){
rlm@23 160
rlm@23 161 dsamples[Points/2-1-i] = dsamples[Points/2-1-i-shift];
rlm@23 162 dimag[Points/2-1-i] = dimag[Points/2-1-i-shift];
rlm@23 163
rlm@23 164
rlm@23 165 dsamples[Points/2+i] = dsamples[Points/2+i+shift];
rlm@23 166 dimag[Points/2+i] = dimag[Points/2+i+shift];
rlm@23 167
rlm@23 168 }
rlm@23 169
rlm@23 170 // fill in the lower points
rlm@23 171 for(int i = 0; i < shift; i++){
rlm@23 172
rlm@23 173 dsamples[i] = 0.0;
rlm@23 174 dimag[i] = 0.0;
rlm@23 175
rlm@23 176 dsamples[Points-1-i] = 0.0;
rlm@23 177 dimag[Points-1-i] = 0.0;
rlm@23 178 }
rlm@23 179
rlm@23 180 DFT(-1,Points,dsamples,dimag);
rlm@23 181
rlm@23 182 short int ifftResult[Points];
rlm@23 183 for(int i = 0; i < Points; i++)
rlm@23 184 ifftResult[i] = (int)dsamples[i];
rlm@23 185
rlm@23 186 // strip off the padding
rlm@23 187 int write = (read < Points/2) ? read+(Points/2) : Points;
rlm@23 188 for(int i = 0; i < write; i++)
rlm@23 189 assert(fwrite(&ifftResult[i],short_sz,1,ifftPcmFile));
rlm@23 190
rlm@23 191 // break out if we're at the end
rlm@23 192 if(read < Points/2){
rlm@23 193 break;
rlm@23 194 }
rlm@23 195 }
rlm@23 196
rlm@23 197 fclose(firPcmFile);
rlm@23 198 fclose(ifftPcmFile);
rlm@23 199
rlm@23 200 // do the windowing
rlm@23 201
rlm@23 202 ifftPcmFile = fopen("ifft.pcm", "r");
rlm@23 203 outputPcmFile = fopen("output.pcm", "w");
rlm@23 204
rlm@23 205 assert(ifftPcmFile);
rlm@23 206 assert(outputPcmFile);
rlm@23 207
rlm@23 208 short int samplesA[Points];
rlm@23 209 short int samplesB[Points];
rlm@23 210
rlm@23 211 bool valid_stream = true;
rlm@23 212 // write out the first half frame
rlm@23 213 for(int i = 0; i < Points/2; i++){
rlm@23 214 valid_stream &= fread(samplesA, short_sz, 1,ifftPcmFile);
rlm@23 215 valid_stream &= fwrite(samplesA, short_sz,1,outputPcmFile);
rlm@23 216 }
rlm@23 217
rlm@23 218 while(valid_stream){
rlm@23 219
rlm@23 220 for(int i = 0; i < Points/2; i++)
rlm@23 221 valid_stream &= fread(&samplesA[i], short_sz, 1,ifftPcmFile);
rlm@23 222
rlm@23 223 for(int i = 0; i < Points; i++)
rlm@23 224 valid_stream &= fread(&samplesB[i], short_sz, 1,ifftPcmFile);
rlm@23 225
rlm@23 226 for(int i = 0; i < Points/2; i++)
rlm@23 227 valid_stream &= fread(&samplesA[(Points/2)+i], short_sz, 1,ifftPcmFile);
rlm@23 228
rlm@23 229
rlm@23 230 // this isn't quite right
rlm@23 231 if(!valid_stream)
rlm@23 232 break;
rlm@23 233
rlm@23 234 for(int i = 0; i < Points; i++){
rlm@23 235 double window = windowing_fn(i);
rlm@23 236 short int rv = (short int)(((double)samplesA[i])*(1.0-window) + window*((double)samplesB[i]));
rlm@23 237 assert(fwrite(&rv,short_sz,1,outputPcmFile));
rlm@23 238 }
rlm@23 239
rlm@23 240 }
rlm@23 241
rlm@23 242 fclose(ifftPcmFile);
rlm@23 243 fclose(outputPcmFile);
rlm@23 244
rlm@23 245 generate_wav("output.pcm", inputWavFileName, outputWavFileName);
rlm@23 246
rlm@23 247 }
rlm@23 248