view tools/audio_processor_test/cosine_window/checker/checker.cpp @ 42:ced2ebd41347 pygar svn.43

[svn r43] bunch of updates that almost work...
author punk
date Wed, 05 May 2010 01:09:09 -0400
parents 90197e3375e2
children
line wrap: on
line source
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include <assert.h>
5 #include <string.h>
6 #include <math.h>
7 #include "SndfileWavUtil.h"
9 const int Taps = 9;
10 const int Points = 8;
11 const int log_Points = 3;
12 static int newest;
13 static short int x[Taps] = {0,0,0,0,0,0,0,0,0};
14 static short int fir(short int);
16 // foreward declaration. This is defined in DFT.cpp
17 int DFT(int dir,int m,double *x1,double *y1);
20 static int x_to_y(int x, int y)
21 {
22 assert(y>0);
23 int rv = x;
24 while(y-->1)
25 x *= x;
26 return rv;
27 }
29 // function implemented in LinearWindowingFunction.bsv
30 static double windowing_fn(int index)
31 {
32 double pi = 3.14159265;
33 double i = (double)index;
34 return sin(0.5*pi* (sin((i + 0.5)/Points * pi) * sin((i + 0.5)/Points * pi)));
35 }
37 // lifted from the bluespec code
38 static double h[Taps] = {-0.0124,
39 0.0,
40 -0.0133,
41 0.0,
42 0.8181,
43 0.0,
44 -0.0133,
45 0.0,
46 -0.0124};
49 short int
50 fir(short int sample)
51 {
52 x[newest] = sample;
53 double y = 0;
54 for (int k = 0; k < Taps; k++) {
55 short int bits = x[(newest+k)%Taps];
56 double x = (double)bits;
57 y += h[k] * x;
58 }
59 newest = (newest == 0) ? (Taps-1) : (newest-1);
60 short int rv = ((short int)y)&0x0000FFFF;
61 return rv;
62 }
65 int
66 main (int argc, char * argv [])
67 {
68 const char* inputWavFileName;
69 const char* outputWavFileName;
71 FILE *inputPcmFile;
72 FILE *firPcmFile;
73 FILE *ifftPcmFile;
74 FILE *outputPcmFile;
76 short int sample;
77 short int samples[Points];
78 const unsigned int short_sz = sizeof(short int);
80 inputWavFileName = argv[1];
81 outputWavFileName = argv[2];
83 // Convert input wav to pcm
84 generate_pcm(inputWavFileName, "input.pcm");
86 inputPcmFile = fopen("input.pcm", "r");
87 firPcmFile = fopen("fir.pcm", "w");
89 assert(inputPcmFile);
90 assert(firPcmFile);
92 newest = 0;
93 memset(x, 0, sizeof(x));
95 while(fread(&sample, short_sz, 1, inputPcmFile)) {
96 sample = fir(sample);
97 assert(fwrite(&sample,short_sz, 1, firPcmFile));
98 }
100 fclose(inputPcmFile);
101 fclose(firPcmFile);
103 firPcmFile = fopen("fir.pcm", "r");
104 ifftPcmFile = fopen("ifft.pcm", "w");
106 assert(firPcmFile);
107 assert(ifftPcmFile);
109 int read = 0;
111 // read in half a frame
112 for(int i = 0; i < Points/2; i++)
113 read += fread(&samples[(Points/2)+i], short_sz, 1, firPcmFile);
115 if(read > 0)
116 assert(read==Points/2);
118 // we will do an 'Points' point fft, and then undo it
119 while(true) {
121 read = 0;
123 // shift last half frame to head
124 for(int i = 0; i < Points/2; i++)
125 samples[i] = samples[(Points/2)+i];
127 // read in another half frame
128 for(int i = 0; i < Points/2; i++)
129 read += fread(&samples[(Points/2)+i], short_sz, 1, firPcmFile);
131 if(read == 0)
132 break;
134 // pad out the
135 if(read < Points/2){
136 for(int i = read; i < Points/2; i++)
137 samples[Points/2+i] = 0;
138 }
140 double dsamples[Points];
141 double dimag[Points];
143 // this shift is performed in the Bluespec
144 for(int i = 0; i < Points; i++){
145 dsamples[i] = (double)(samples[i]>>log_Points);
146 dimag[i] = 0.0;
147 }
150 //for(int i = 0; i < Points; i++)
151 // printf("%d ", (int)dsamples[i]);
152 //printf("\n");
154 DFT(1,Points,dsamples,dimag);
156 // shift up by 'shift' (in)
157 int shift = 2;
159 for(int i = 0; i < Points/2; i++){
161 dsamples[Points/2-1-i] = dsamples[Points/2-1-i-shift];
162 dimag[Points/2-1-i] = dimag[Points/2-1-i-shift];
165 dsamples[Points/2+i] = dsamples[Points/2+i+shift];
166 dimag[Points/2+i] = dimag[Points/2+i+shift];
168 }
170 // fill in the lower points
171 for(int i = 0; i < shift; i++){
173 dsamples[i] = 0.0;
174 dimag[i] = 0.0;
176 dsamples[Points-1-i] = 0.0;
177 dimag[Points-1-i] = 0.0;
178 }
180 DFT(-1,Points,dsamples,dimag);
182 short int ifftResult[Points];
183 for(int i = 0; i < Points; i++)
184 ifftResult[i] = (int)dsamples[i];
186 // strip off the padding
187 int write = (read < Points/2) ? read+(Points/2) : Points;
188 for(int i = 0; i < write; i++)
189 assert(fwrite(&ifftResult[i],short_sz,1,ifftPcmFile));
191 // break out if we're at the end
192 if(read < Points/2){
193 break;
194 }
195 }
197 fclose(firPcmFile);
198 fclose(ifftPcmFile);
200 // do the windowing
202 ifftPcmFile = fopen("ifft.pcm", "r");
203 outputPcmFile = fopen("output.pcm", "w");
205 assert(ifftPcmFile);
206 assert(outputPcmFile);
208 short int samplesA[Points];
209 short int samplesB[Points];
211 bool valid_stream = true;
212 // write out the first half frame
213 for(int i = 0; i < Points/2; i++){
214 valid_stream &= fread(samplesA, short_sz, 1,ifftPcmFile);
215 valid_stream &= fwrite(samplesA, short_sz,1,outputPcmFile);
216 }
218 while(valid_stream){
220 for(int i = 0; i < Points/2; i++)
221 valid_stream &= fread(&samplesA[i], short_sz, 1,ifftPcmFile);
223 for(int i = 0; i < Points; i++)
224 valid_stream &= fread(&samplesB[i], short_sz, 1,ifftPcmFile);
226 for(int i = 0; i < Points/2; i++)
227 valid_stream &= fread(&samplesA[(Points/2)+i], short_sz, 1,ifftPcmFile);
230 // this isn't quite right
231 if(!valid_stream)
232 break;
234 for(int i = 0; i < Points; i++){
235 double window = windowing_fn(i);
236 short int rv = (short int)(((double)samplesA[i])*(1.0-window) + window*((double)samplesB[i]));
237 assert(fwrite(&rv,short_sz,1,outputPcmFile));
238 }
240 }
242 fclose(ifftPcmFile);
243 fclose(outputPcmFile);
245 generate_wav("output.pcm", inputWavFileName, outputWavFileName);
247 }