Mercurial > pygar

comparison tools/audio_processor_test/transform/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 "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
19 static int x_to_y(int x, int y)
20 {
21 assert(y>0);
22 int rv = x;
23 while(y-->1)
24 x *= x;
25 return rv;
26 }
27
28 // function implemented in LinearWindowingFunction.bsv
29 static double windowing_fn(int index)
30 {
31 double i = (double)index;
32 double divisor = Points/2; //(double) (x_to_y(2,log_Points-1));
33 double rv = (index < Points/2) ? (i / divisor) : ((((double)Points)-i) / divisor);
34 return rv;
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 DFT(1,Points,dsamples,dimag);
151 DFT(-1,Points,dsamples,dimag);
152
153 short int ifftResult[Points];
154 for(int i = 0; i < Points; i++)
155 ifftResult[i] = (int)dsamples[i]; //)>>log_Points;
156
157 // strip off the padding
158 int write = (read < Points/2) ? read+(Points/2) : Points;
159 for(int i = 0; i < write; i++)
160 assert(fwrite(&ifftResult[i],short_sz,1,ifftPcmFile));
161
162 // break out if we're at the end
163 if(read < Points/2){
164 break;
165 }
166 }
167
168 fclose(firPcmFile);
169 fclose(ifftPcmFile);
170
171 // do the windowing
172
173 ifftPcmFile = fopen("ifft.pcm", "r");
174 outputPcmFile = fopen("output.pcm", "w");
175
176 assert(ifftPcmFile);
177 assert(outputPcmFile);
178
179 short int samplesA[Points];
180 short int samplesB[Points];
181
182 bool valid_stream = true;
183 // write out the first half frame
184 for(int i = 0; i < Points/2; i++){
185 valid_stream &= fread(samplesA, short_sz, 1,ifftPcmFile);
186 valid_stream &= fwrite(samplesA, short_sz,1,outputPcmFile);
187 }
188
189 while(valid_stream){
190
191 for(int i = 0; i < Points/2; i++)
192 valid_stream &= fread(&samplesA[i], short_sz, 1,ifftPcmFile);
193
194 for(int i = 0; i < Points; i++)
195 valid_stream &= fread(&samplesB[i], short_sz, 1,ifftPcmFile);
196
197 for(int i = 0; i < Points/2; i++)
198 valid_stream &= fread(&samplesA[(Points/2)+i], short_sz, 1,ifftPcmFile);
199
200
201 // this isn't quite right
202 if(!valid_stream)
203 break;
204
205 for(int i = 0; i < Points; i++){
206 double window = windowing_fn(i);
207 short int rv = samplesA[i]*(1-window) + window*samplesB[i];
208 assert(fwrite(&rv,short_sz,1,outputPcmFile));
209 }
210
211 }
212
213 fclose(ifftPcmFile);
214 fclose(outputPcmFile);
215
216 generate_wav("output.pcm", inputWavFileName, outputWavFileName);
217
218 }
219