Mercurial > pygar

comparison tools/audio_processor_test/cosine_window/checker/checker.cpp @ 23:90197e3375e2 pygar svn.24

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[svn r24] added testing, but something is wrong with our c++ file.
author rlm
date Wed, 28 Apr 2010 08:19:09 -0400
parents
children
comparison
equal deleted inserted replaced
22:0cfbb1e2de22 23:90197e3375e2
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"
8
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);
15
16 // foreward declaration. This is defined in DFT.cpp
17 int DFT(int dir,int m,double *x1,double *y1);
18
19
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 }
28
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 }
36
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};
47
48
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 }
63
64
65 int
66 main (int argc, char * argv [])
67 {
68 const char* inputWavFileName;
69 const char* outputWavFileName;
70
71 FILE *inputPcmFile;
72 FILE *firPcmFile;
73 FILE *ifftPcmFile;
74 FILE *outputPcmFile;
75
76 short int sample;
77 short int samples[Points];
78 const unsigned int short_sz = sizeof(short int);
79
80 inputWavFileName = argv[1];
81 outputWavFileName = argv[2];
82
83 // Convert input wav to pcm
84 generate_pcm(inputWavFileName, "input.pcm");
85
86 inputPcmFile = fopen("input.pcm", "r");
87 firPcmFile = fopen("fir.pcm", "w");
88
89 assert(inputPcmFile);
90 assert(firPcmFile);
91
92 newest = 0;
93 memset(x, 0, sizeof(x));
94
95 while(fread(&sample, short_sz, 1, inputPcmFile)) {
96 sample = fir(sample);
97 assert(fwrite(&sample,short_sz, 1, firPcmFile));
98 }
99
100 fclose(inputPcmFile);
101 fclose(firPcmFile);
102
103 firPcmFile = fopen("fir.pcm", "r");
104 ifftPcmFile = fopen("ifft.pcm", "w");
105
106 assert(firPcmFile);
107 assert(ifftPcmFile);
108
109 int read = 0;
110
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);
114
115 if(read > 0)
116 assert(read==Points/2);
117
118 // we will do an 'Points' point fft, and then undo it
119 while(true) {
120
121 read = 0;
122
123 // shift last half frame to head
124 for(int i = 0; i < Points/2; i++)
125 samples[i] = samples[(Points/2)+i];
126
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);
130
131 if(read == 0)
132 break;
133
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 }
139
140 double dsamples[Points];
141 double dimag[Points];
142
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 }
148
149
150 //for(int i = 0; i < Points; i++)
151 // printf("%d ", (int)dsamples[i]);
152 //printf("\n");
153
154 DFT(1,Points,dsamples,dimag);
155
156 // shift up by 'shift' (in)
157 int shift = 2;
158
159 for(int i = 0; i < Points/2; i++){
160
161 dsamples[Points/2-1-i] = dsamples[Points/2-1-i-shift];
162 dimag[Points/2-1-i] = dimag[Points/2-1-i-shift];
163
164
165 dsamples[Points/2+i] = dsamples[Points/2+i+shift];
166 dimag[Points/2+i] = dimag[Points/2+i+shift];
167
168 }
169
170 // fill in the lower points
171 for(int i = 0; i < shift; i++){
172
173 dsamples[i] = 0.0;
174 dimag[i] = 0.0;
175
176 dsamples[Points-1-i] = 0.0;
177 dimag[Points-1-i] = 0.0;
178 }
179
180 DFT(-1,Points,dsamples,dimag);
181
182 short int ifftResult[Points];
183 for(int i = 0; i < Points; i++)
184 ifftResult[i] = (int)dsamples[i];
185
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));
190
191 // break out if we're at the end
192 if(read < Points/2){
193 break;
194 }
195 }
196
197 fclose(firPcmFile);
198 fclose(ifftPcmFile);
199
200 // do the windowing
201
202 ifftPcmFile = fopen("ifft.pcm", "r");
203 outputPcmFile = fopen("output.pcm", "w");
204
205 assert(ifftPcmFile);
206 assert(outputPcmFile);
207
208 short int samplesA[Points];
209 short int samplesB[Points];
210
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 }
217
218 while(valid_stream){
219
220 for(int i = 0; i < Points/2; i++)
221 valid_stream &= fread(&samplesA[i], short_sz, 1,ifftPcmFile);
222
223 for(int i = 0; i < Points; i++)
224 valid_stream &= fread(&samplesB[i], short_sz, 1,ifftPcmFile);
225
226 for(int i = 0; i < Points/2; i++)
227 valid_stream &= fread(&samplesA[(Points/2)+i], short_sz, 1,ifftPcmFile);
228
229
230 // this isn't quite right
231 if(!valid_stream)
232 break;
233
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 }
239
240 }
241
242 fclose(ifftPcmFile);
243 fclose(outputPcmFile);
244
245 generate_wav("output.pcm", inputWavFileName, outputWavFileName);
246
247 }
248