15 #define MUL16(a,b) ((a) * (b))
17 #define CMAC(pre, pim, are, aim, bre, bim) \
19 pre += (MUL16(are, bre) - MUL16(aim, bim));\
20 pim += (MUL16(are, bim) + MUL16(bre, aim));\
25 void fft_ref_init(int nbits
, int inverse
)
31 exptab
= av_malloc((n
/ 2) * sizeof(FFTComplex
));
33 for(i
=0;i
<(n
/2);i
++) {
34 alpha
= 2 * M_PI
* (float)i
/ (float)n
;
44 void fft_ref(FFTComplex
*tabr
, FFTComplex
*tab
, int nbits
)
47 float tmp_re
, tmp_im
, s
, c
;
57 k
= (i
* j
) & (n
- 1);
59 c
= -exptab
[k
- n2
].re
;
60 s
= -exptab
[k
- n2
].im
;
65 CMAC(tmp_re
, tmp_im
, c
, s
, q
->re
, q
->im
);
73 void imdct_ref(float *out
, float *in
, int n
)
81 a
= (2 * i
+ 1 + (n
/ 2)) * (2 * k
+ 1);
82 f
= cos(M_PI
* a
/ (double)(2 * n
));
89 /* NOTE: no normalisation by 1 / N is done */
90 void mdct_ref(float *output
, float *input
, int n
)
99 a
= (2*M_PI
*(2*i
+1+n
/2)*(2*k
+1) / (4 * n
));
100 s
+= input
[i
] * cos(a
);
109 return (float)((random() & 0xffff) - 32768) / 32768.0;
112 int64_t gettime(void)
115 gettimeofday(&tv
,NULL
);
116 return (int64_t)tv
.tv_sec
* 1000000 + tv
.tv_usec
;
119 void check_diff(float *tab1
, float *tab2
, int n
)
124 if (fabsf(tab1
[i
] - tab2
[i
]) >= 1e-3) {
125 av_log(NULL
, AV_LOG_ERROR
, "ERROR %d: %f %f\n",
126 i
, tab1
[i
], tab2
[i
]);
134 av_log(NULL
, AV_LOG_INFO
,"usage: fft-test [-h] [-s] [-i] [-n b]\n"
135 "-h print this help\n"
138 "-i inverse transform test\n"
139 "-n b set the transform size to 2^b\n"
146 int main(int argc
, char **argv
)
148 FFTComplex
*tab
, *tab1
, *tab_ref
;
149 FFTSample
*tabtmp
, *tab2
;
154 FFTContext s1
, *s
= &s1
;
155 MDCTContext m1
, *m
= &m1
;
156 int fft_nbits
, fft_size
;
161 c
= getopt(argc
, argv
, "hsimn:");
178 fft_nbits
= atoi(optarg
);
183 fft_size
= 1 << fft_nbits
;
184 tab
= av_malloc(fft_size
* sizeof(FFTComplex
));
185 tab1
= av_malloc(fft_size
* sizeof(FFTComplex
));
186 tab_ref
= av_malloc(fft_size
* sizeof(FFTComplex
));
187 tabtmp
= av_malloc(fft_size
/ 2 * sizeof(FFTSample
));
188 tab2
= av_malloc(fft_size
* sizeof(FFTSample
));
192 av_log(NULL
, AV_LOG_INFO
,"IMDCT");
194 av_log(NULL
, AV_LOG_INFO
,"MDCT");
195 ff_mdct_init(m
, fft_nbits
, do_inverse
);
198 av_log(NULL
, AV_LOG_INFO
,"IFFT");
200 av_log(NULL
, AV_LOG_INFO
,"FFT");
201 ff_fft_init(s
, fft_nbits
, do_inverse
);
202 fft_ref_init(fft_nbits
, do_inverse
);
204 av_log(NULL
, AV_LOG_INFO
," %d test\n", fft_size
);
206 /* generate random data */
208 for(i
=0;i
<fft_size
;i
++) {
209 tab1
[i
].re
= frandom();
210 tab1
[i
].im
= frandom();
213 /* checking result */
214 av_log(NULL
, AV_LOG_INFO
,"Checking...\n");
218 imdct_ref((float *)tab_ref
, (float *)tab1
, fft_size
);
219 ff_imdct_calc(m
, tab2
, (float *)tab1
, tabtmp
);
220 check_diff((float *)tab_ref
, tab2
, fft_size
);
222 mdct_ref((float *)tab_ref
, (float *)tab1
, fft_size
);
224 ff_mdct_calc(m
, tab2
, (float *)tab1
, tabtmp
);
226 check_diff((float *)tab_ref
, tab2
, fft_size
/ 2);
229 memcpy(tab
, tab1
, fft_size
* sizeof(FFTComplex
));
230 ff_fft_permute(s
, tab
);
233 fft_ref(tab_ref
, tab1
, fft_nbits
);
234 check_diff((float *)tab_ref
, (float *)tab
, fft_size
* 2);
237 /* do a speed test */
240 int64_t time_start
, duration
;
243 av_log(NULL
, AV_LOG_INFO
,"Speed test...\n");
244 /* we measure during about 1 seconds */
247 time_start
= gettime();
248 for(it
=0;it
<nb_its
;it
++) {
251 ff_imdct_calc(m
, (float *)tab
, (float *)tab1
, tabtmp
);
253 ff_mdct_calc(m
, (float *)tab
, (float *)tab1
, tabtmp
);
256 memcpy(tab
, tab1
, fft_size
* sizeof(FFTComplex
));
260 duration
= gettime() - time_start
;
261 if (duration
>= 1000000)
265 av_log(NULL
, AV_LOG_INFO
,"time: %0.1f us/transform [total time=%0.2f s its=%d]\n",
266 (double)duration
/ nb_its
,
267 (double)duration
/ 1000000.0,