44d19a651743754ba5114942314d21db1cfea893
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <math.h>
5 #include <sys/time.h>
7 #include <resample.h>
9 #define AMP 16000
10 #define I_RATE 48000
11 #define O_RATE 44100
12 //#define O_RATE 24000
14 //#define test_func(x) 1
15 //#define test_func(x) sin(2*M_PI*(x)*10)
16 //#define test_func(x) sin(2*M_PI*(x)*(x)*1000)
17 #define test_func(x) sin(2*M_PI*(x)*(x)*12000)
19 short i_buf[I_RATE*2*2];
20 short o_buf[O_RATE*2*2];
22 static int i_offset;
23 static int o_offset;
25 FILE *out;
27 void test_res1(void);
28 void test_res2(void);
29 void test_res3(void);
30 void test_res4(void);
31 void test_res5(void);
32 void test_res6(void);
33 void test_res7(void);
35 int main(int argc,char *argv[])
36 {
37 out = fopen("out","w");
39 test_res7();
41 return 0;
42 }
44 void *get_buffer(void *priv, unsigned int size)
45 {
46 void *ret;
47 ret = ((void *)o_buf) + o_offset;
48 o_offset += size;
49 return ret;
50 }
52 struct timeval start_time;
53 void start_timer(void)
54 {
55 gettimeofday(&start_time,NULL);
56 //printf("start %ld.%06ld\n",start_time.tv_sec,start_time.tv_usec);
57 }
59 void end_timer(void)
60 {
61 struct timeval end_time;
62 double diff;
64 gettimeofday(&end_time,NULL);
65 //printf("end %ld.%06ld\n",end_time.tv_sec,end_time.tv_usec);
66 diff = (end_time.tv_sec - start_time.tv_sec) +
67 1e-6*(end_time.tv_usec - start_time.tv_usec);
69 printf("time %g\n",diff);
71 }
73 void test_res1(void)
74 {
75 resample_t *r;
76 int i;
77 double sum10k,sum22k;
78 double f;
79 int n10k,n22k;
80 double x;
82 for(i=0;i<I_RATE;i++){
83 i_buf[i*2+0] = rint(AMP * test_func((double)i/I_RATE));
84 //i_buf[i*2+1] = rint(AMP * test_func((double)i/I_RATE));
85 i_buf[i*2+1] = (i<1000)?AMP:0;
86 }
88 r = malloc(sizeof(resample_t));
89 memset(r,0,sizeof(resample_t));
91 r->i_rate = I_RATE;
92 r->o_rate = O_RATE;
93 //r->method = RESAMPLE_SINC_SLOW;
94 r->method = RESAMPLE_SINC;
95 r->channels = 2;
96 //r->verbose = 1;
97 r->filter_length = 64;
98 r->get_buffer = get_buffer;
100 resample_init(r);
102 start_timer();
103 #define blocked
104 #ifdef blocked
105 for(i=0;i+256<I_RATE;i+=256){
106 resample_scale(r,i_buf+i*2,256*2*2);
107 }
108 if(I_RATE-i){
109 resample_scale(r,i_buf+i*2,(I_RATE-i)*2*2);
110 }
111 #else
112 resample_scale(r,i_buf,I_RATE*2*2);
113 #endif
114 end_timer();
116 for(i=0;i<O_RATE;i++){
117 f = AMP*test_func((double)i/O_RATE);
118 //f = rint(AMP*test_func((double)i/O_RATE));
119 fprintf(out,"%d %d %d %g %g\n",i,
120 o_buf[2*i+0],o_buf[2*i+1],
121 f,o_buf[2*i+0]-f);
122 }
124 sum10k=0;
125 sum22k=0;
126 n10k=0;
127 n22k=0;
128 for(i=0;i<O_RATE;i++){
129 f = AMP*test_func((double)i/O_RATE);
130 //f = rint(AMP*test_func((double)i/O_RATE));
131 x = o_buf[2*i+0]-f;
132 if(((0.5*i)/O_RATE*I_RATE)<10000){
133 sum10k += x*x;
134 n10k++;
135 }
136 if(((0.5*i)/O_RATE*I_RATE)<22050){
137 sum22k += x*x;
138 n22k++;
139 }
140 }
141 printf("average error 10k=%g 22k=%g\n",
142 sqrt(sum10k/n10k),
143 sqrt(sum22k/n22k));
144 }
147 void test_res2(void)
148 {
149 functable_t *t;
150 int i;
151 double x;
152 double f1,f2;
154 t = malloc(sizeof(*t));
155 memset(t,0,sizeof(*t));
157 t->start = -50.0;
158 t->offset = 1;
159 t->len = 100;
161 t->func_x = functable_sinc;
162 t->func_dx = functable_dsinc;
164 functable_init(t);
166 for(i=0;i<1000;i++){
167 x = -50.0 + 0.1 * i;
168 f1 = functable_sinc(NULL,x);
169 f2 = functable_eval(t,x);
170 fprintf(out,"%d %g %g %g\n",i,f1,f2,f1-f2);
171 }
172 }
174 void test_res3(void)
175 {
176 functable_t *t;
177 int i;
178 double x;
179 double f1,f2;
180 int n = 1;
182 t = malloc(sizeof(*t));
183 memset(t,0,sizeof(*t));
185 t->start = -50.0;
186 t->offset = 1.0 / n;
187 t->len = 100 * n;
189 t->func_x = functable_sinc;
190 t->func_dx = functable_dsinc;
192 t->func2_x = functable_window_std;
193 t->func2_dx = functable_window_dstd;
195 t->scale = 1.0;
196 t->scale2 = 1.0 / (M_PI * 16);
198 functable_init(t);
200 for(i=0;i<1000 * n;i++){
201 x = -50.0 + 0.1/n * i;
202 f1 = functable_sinc(NULL,t->scale * x) *
203 functable_window_std(NULL,t->scale2 * x);
204 f2 = functable_eval(t,x);
205 fprintf(out,"%d %g %g %g\n",i,f1,f2,f2-f1);
206 }
207 }
209 double sinc_poly(double x)
210 {
211 #define INV3FAC 1.66666666666666666e-1
212 #define INV5FAC 8.33333333333333333e-3
213 #define INV7FAC 1.984126984e-4
214 #define INV9FAC 2.755731922e-6
215 #define INV11FAC 2.505210839e-8
216 double x2 = x * x;
218 return 1
219 - x2 * INV3FAC
220 + x2 * x2 * INV5FAC
221 - x2 * x2 * x2 * INV7FAC;
222 //+ x2 * x2 * x2 * x2 * INV9FAC
223 //- x2 * x2 * x2 * x2 * x2 * INV11FAC;
224 }
226 void test_res4(void)
227 {
228 int i;
229 double x,f1,f2;
231 for(i=1;i<100;i++){
232 x = 0.01 * i;
233 f1 = 1 - sin(x)/x;
234 f2 = 1 - sinc_poly(x);
236 fprintf(out,"%g %.20g %.20g %.20g\n",x,f1,f2,f2-f1);
237 }
238 }
241 void test_res5(void)
242 {
243 int i;
244 double sum;
246 start_timer();
247 sum = 0;
248 for(i=0;i<I_RATE;i++){
249 sum += i_buf[i*2];
250 }
251 end_timer();
252 i_buf[0] = sum;
253 }
256 void short_to_double(double *d,short *x) { *d = *x; }
257 void short_to_float(float *f,short *x) { *f = *x; }
258 void float_to_double(double *f,float *x) { *f = *x; }
259 void double_to_short(short *f,double *x) { *f = *x; }
261 double res6_tmp[1000];
263 void test_res6(void)
264 {
265 int i;
267 for(i=0;i<I_RATE;i++){
268 i_buf[i] = rint(AMP * test_func((double)i/I_RATE));
269 }
271 conv_double_short_ref(res6_tmp,i_buf,1000);
272 for(i=0;i<1000;i++){
273 res6_tmp[i] *= 3.0;
274 }
275 conv_short_double_ppcasm(o_buf,res6_tmp,1000);
277 for(i=0;i<1000;i++){
278 fprintf(out,"%d %d %g %d\n",i,i_buf[i],res6_tmp[i],o_buf[i]);
279 }
280 }
282 void test_res7(void)
283 {
284 resample_t *r;
285 int i;
286 double sum10k,sum22k;
287 double f;
288 int n10k,n22k;
289 double x;
291 for(i=0;i<I_RATE;i++){
292 i_buf[i] = rint(AMP * test_func((double)i/I_RATE));
293 }
295 r = malloc(sizeof(resample_t));
296 memset(r,0,sizeof(resample_t));
298 r->i_rate = I_RATE;
299 r->o_rate = O_RATE;
300 //r->method = RESAMPLE_SINC_SLOW;
301 r->method = RESAMPLE_SINC;
302 r->channels = 1;
303 //r->verbose = 1;
304 r->filter_length = 64;
305 r->get_buffer = get_buffer;
307 resample_init(r);
309 start_timer();
310 #define blocked
311 #ifdef blocked
312 for(i=0;i+256<I_RATE;i+=256){
313 resample_scale(r,i_buf+i,256*2);
314 }
315 if(I_RATE-i){
316 resample_scale(r,i_buf+i,(I_RATE-i)*2);
317 }
318 #else
319 resample_scale(r,i_buf,I_RATE*2);
320 #endif
321 end_timer();
323 for(i=0;i<O_RATE;i++){
324 f = AMP*test_func((double)i/O_RATE);
325 //f = rint(AMP*test_func((double)i/O_RATE));
326 fprintf(out,"%d %d %d %g %g\n",i,
327 o_buf[i],0,
328 f,o_buf[i]-f);
329 }
331 sum10k=0;
332 sum22k=0;
333 n10k=0;
334 n22k=0;
335 for(i=0;i<O_RATE;i++){
336 f = AMP*test_func((double)i/O_RATE);
337 //f = rint(AMP*test_func((double)i/O_RATE));
338 x = o_buf[i]-f;
339 if(((0.5*i)/O_RATE*I_RATE)<10000){
340 sum10k += x*x;
341 n10k++;
342 }
343 if(((0.5*i)/O_RATE*I_RATE)<22050){
344 sum22k += x*x;
345 n22k++;
346 }
347 }
348 printf("average error 10k=%g 22k=%g\n",
349 sqrt(sum10k/n10k),
350 sqrt(sum22k/n22k));
351 }