## ffmpeg / libavcodec / fft-test.c @ d4c5d2ad

History | View | Annotate | Download (6.28 KB)

1 | 983e3246 | Michael Niedermayer | ```
/**
``` |
---|---|---|---|

2 | ```
* @file fft-test.c
``` |
||

3 | ```
* FFT and MDCT tests.
``` |
||

4 | ```
*/
``` |
||

5 | |||

6 | bb6f5690 | Fabrice Bellard | #include "dsputil.h" |

7 | #include <math.h> |
||

8 | e366e679 | Fabrice Bellard | #include <unistd.h> |

9 | bb6f5690 | Fabrice Bellard | #include <sys/time.h> |

10 | |||

11 | ```
int mm_flags;
``` |
||

12 | |||

13 | ```
/* reference fft */
``` |
||

14 | |||

15 | ```
#define MUL16(a,b) ((a) * (b))
``` |
||

16 | |||

17 | ```
#define CMAC(pre, pim, are, aim, bre, bim) \
``` |
||

18 | {\ |
||

19 | pre += (MUL16(are, bre) - MUL16(aim, bim));\ |
||

20 | pim += (MUL16(are, bim) + MUL16(bre, aim));\ |
||

21 | } |
||

22 | |||

23 | FFTComplex *exptab; |
||

24 | |||

25 | void fft_ref_init(int nbits, int inverse) |
||

26 | { |
||

27 | ```
int n, i;
``` |
||

28 | ```
float c1, s1, alpha;
``` |
||

29 | |||

30 | ```
n = 1 << nbits;
``` |
||

31 | exptab = av_malloc((n / 2) * sizeof(FFTComplex)); |
||

32 | |||

33 | for(i=0;i<(n/2);i++) { |
||

34 | alpha = 2 * M_PI * (float)i / (float)n; |
||

35 | c1 = cos(alpha); |
||

36 | s1 = sin(alpha); |
||

37 | ```
if (!inverse)
``` |
||

38 | s1 = -s1; |
||

39 | exptab[i].re = c1; |
||

40 | exptab[i].im = s1; |
||

41 | } |
||

42 | } |
||

43 | |||

44 | void fft_ref(FFTComplex *tabr, FFTComplex *tab, int nbits) |
||

45 | { |
||

46 | ```
int n, i, j, k, n2;
``` |
||

47 | ```
float tmp_re, tmp_im, s, c;
``` |
||

48 | FFTComplex *q; |
||

49 | |||

50 | ```
n = 1 << nbits;
``` |
||

51 | ```
n2 = n >> 1;
``` |
||

52 | for(i=0;i<n;i++) { |
||

53 | ```
tmp_re = 0;
``` |
||

54 | ```
tmp_im = 0;
``` |
||

55 | q = tab; |
||

56 | for(j=0;j<n;j++) { |
||

57 | ```
k = (i * j) & (n - 1);
``` |
||

58 | ```
if (k >= n2) {
``` |
||

59 | c = -exptab[k - n2].re; |
||

60 | s = -exptab[k - n2].im; |
||

61 | ```
} else {
``` |
||

62 | c = exptab[k].re; |
||

63 | s = exptab[k].im; |
||

64 | } |
||

65 | CMAC(tmp_re, tmp_im, c, s, q->re, q->im); |
||

66 | q++; |
||

67 | } |
||

68 | tabr[i].re = tmp_re; |
||

69 | tabr[i].im = tmp_im; |
||

70 | } |
||

71 | } |
||

72 | |||

73 | void imdct_ref(float *out, float *in, int n) |
||

74 | { |
||

75 | ```
int k, i, a;
``` |
||

76 | ```
float sum, f;
``` |
||

77 | |||

78 | for(i=0;i<n;i++) { |
||

79 | ```
sum = 0;
``` |
||

80 | for(k=0;k<n/2;k++) { |
||

81 | a = (2 * i + 1 + (n / 2)) * (2 * k + 1); |
||

82 | f = cos(M_PI * a / (double)(2 * n)); |
||

83 | sum += f * in[k]; |
||

84 | } |
||

85 | out[i] = -sum; |
||

86 | } |
||

87 | } |
||

88 | |||

89 | ```
/* NOTE: no normalisation by 1 / N is done */
``` |
||

90 | void mdct_ref(float *output, float *input, int n) |
||

91 | { |
||

92 | ```
int k, i;
``` |
||

93 | ```
float a, s;
``` |
||

94 | |||

95 | ```
/* do it by hand */
``` |
||

96 | for(k=0;k<n/2;k++) { |
||

97 | ```
s = 0;
``` |
||

98 | for(i=0;i<n;i++) { |
||

99 | a = (2*M_PI*(2*i+1+n/2)*(2*k+1) / (4 * n)); |
||

100 | s += input[i] * cos(a); |
||

101 | } |
||

102 | output[k] = s; |
||

103 | } |
||

104 | } |
||

105 | |||

106 | |||

107 | float frandom(void) |
||

108 | { |
||

109 | return (float)((random() & 0xffff) - 32768) / 32768.0; |
||

110 | } |
||

111 | |||

112 | 0c1a9eda | Zdenek Kabelac | ```
int64_t gettime(void)
``` |

113 | bb6f5690 | Fabrice Bellard | { |

114 | ```
struct timeval tv;
``` |
||

115 | ```
gettimeofday(&tv,NULL);
``` |
||

116 | 0c1a9eda | Zdenek Kabelac | return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec; |

117 | bb6f5690 | Fabrice Bellard | } |

118 | |||

119 | void check_diff(float *tab1, float *tab2, int n) |
||

120 | { |
||

121 | ```
int i;
``` |
||

122 | |||

123 | for(i=0;i<n;i++) { |
||

124 | if (fabsf(tab1[i] - tab2[i]) >= 1e-3) { |
||

125 | ```
printf("ERROR %d: %f %f\n",
``` |
||

126 | i, tab1[i], tab2[i]); |
||

127 | } |
||

128 | } |
||

129 | } |
||

130 | |||

131 | |||

132 | void help(void) |
||

133 | { |
||

134 | ```
printf("usage: fft-test [-h] [-s] [-i] [-n b]\n"
``` |
||

135 | ```
"-h print this help\n"
``` |
||

136 | ```
"-s speed test\n"
``` |
||

137 | ```
"-m (I)MDCT test\n"
``` |
||

138 | ```
"-i inverse transform test\n"
``` |
||

139 | ```
"-n b set the transform size to 2^b\n"
``` |
||

140 | ); |
||

141 | ```
exit(1);
``` |
||

142 | } |
||

143 | |||

144 | |||

145 | |||

146 | int main(int argc, char **argv) |
||

147 | { |
||

148 | FFTComplex *tab, *tab1, *tab_ref; |
||

149 | FFTSample *tabtmp, *tab2; |
||

150 | ```
int it, i, c;
``` |
||

151 | int do_speed = 0; |
||

152 | int do_mdct = 0; |
||

153 | int do_inverse = 0; |
||

154 | FFTContext s1, *s = &s1; |
||

155 | MDCTContext m1, *m = &m1; |
||

156 | ```
int fft_nbits, fft_size;
``` |
||

157 | |||

158 | ```
mm_flags = 0;
``` |
||

159 | ```
fft_nbits = 9;
``` |
||

160 | ```
for(;;) {
``` |
||

161 | ```
c = getopt(argc, argv, "hsimn:");
``` |
||

162 | if (c == -1) |
||

163 | ```
break;
``` |
||

164 | ```
switch(c) {
``` |
||

165 | case 'h': |
||

166 | help(); |
||

167 | ```
break;
``` |
||

168 | case 's': |
||

169 | ```
do_speed = 1;
``` |
||

170 | ```
break;
``` |
||

171 | case 'i': |
||

172 | ```
do_inverse = 1;
``` |
||

173 | ```
break;
``` |
||

174 | case 'm': |
||

175 | ```
do_mdct = 1;
``` |
||

176 | ```
break;
``` |
||

177 | case 'n': |
||

178 | fft_nbits = atoi(optarg); |
||

179 | ```
break;
``` |
||

180 | } |
||

181 | } |
||

182 | |||

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));
``` |
||

189 | |||

190 | ```
if (do_mdct) {
``` |
||

191 | ```
if (do_inverse)
``` |
||

192 | ```
printf("IMDCT");
``` |
||

193 | ```
else
``` |
||

194 | ```
printf("MDCT");
``` |
||

195 | ea0f841a | Fabrice Bellard | ff_mdct_init(m, fft_nbits, do_inverse); |

196 | bb6f5690 | Fabrice Bellard | ```
} else {
``` |

197 | ```
if (do_inverse)
``` |
||

198 | ```
printf("IFFT");
``` |
||

199 | ```
else
``` |
||

200 | ```
printf("FFT");
``` |
||

201 | 68951ecf | Gildas Bazin | ff_fft_init(s, fft_nbits, do_inverse); |

202 | bb6f5690 | Fabrice Bellard | fft_ref_init(fft_nbits, do_inverse); |

203 | } |
||

204 | ```
printf(" %d test\n", fft_size);
``` |
||

205 | |||

206 | ```
/* generate random data */
``` |
||

207 | |||

208 | for(i=0;i<fft_size;i++) { |
||

209 | tab1[i].re = frandom(); |
||

210 | tab1[i].im = frandom(); |
||

211 | } |
||

212 | |||

213 | ```
/* checking result */
``` |
||

214 | ```
printf("Checking...\n");
``` |
||

215 | |||

216 | ```
if (do_mdct) {
``` |
||

217 | ```
if (do_inverse) {
``` |
||

218 | imdct_ref((float *)tab_ref, (float *)tab1, fft_size); |
||

219 | ea0f841a | Fabrice Bellard | ```
ff_imdct_calc(m, tab2, (float *)tab1, tabtmp);
``` |

220 | bb6f5690 | Fabrice Bellard | ```
check_diff((float *)tab_ref, tab2, fft_size);
``` |

221 | ```
} else {
``` |
||

222 | mdct_ref((float *)tab_ref, (float *)tab1, fft_size); |
||

223 | |||

224 | ea0f841a | Fabrice Bellard | ```
ff_mdct_calc(m, tab2, (float *)tab1, tabtmp);
``` |

225 | bb6f5690 | Fabrice Bellard | |

226 | check_diff((float *)tab_ref, tab2, fft_size / 2); |
||

227 | } |
||

228 | ```
} else {
``` |
||

229 | ```
memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
``` |
||

230 | 68951ecf | Gildas Bazin | ff_fft_permute(s, tab); |

231 | ff_fft_calc(s, tab); |
||

232 | bb6f5690 | Fabrice Bellard | |

233 | fft_ref(tab_ref, tab1, fft_nbits); |
||

234 | check_diff((float *)tab_ref, (float *)tab, fft_size * 2); |
||

235 | } |
||

236 | |||

237 | ```
/* do a speed test */
``` |
||

238 | |||

239 | ```
if (do_speed) {
``` |
||

240 | 0c1a9eda | Zdenek Kabelac | int64_t time_start, duration; |

241 | bb6f5690 | Fabrice Bellard | ```
int nb_its;
``` |

242 | |||

243 | ```
printf("Speed test...\n");
``` |
||

244 | ```
/* we measure during about 1 seconds */
``` |
||

245 | ```
nb_its = 1;
``` |
||

246 | ```
for(;;) {
``` |
||

247 | time_start = gettime(); |
||

248 | for(it=0;it<nb_its;it++) { |
||

249 | ```
if (do_mdct) {
``` |
||

250 | ```
if (do_inverse) {
``` |
||

251 | ea0f841a | Fabrice Bellard | ff_imdct_calc(m, (float *)tab, (float *)tab1, tabtmp); |

252 | bb6f5690 | Fabrice Bellard | ```
} else {
``` |

253 | ea0f841a | Fabrice Bellard | ff_mdct_calc(m, (float *)tab, (float *)tab1, tabtmp); |

254 | bb6f5690 | Fabrice Bellard | } |

255 | ```
} else {
``` |
||

256 | ```
memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
``` |
||

257 | 68951ecf | Gildas Bazin | ff_fft_calc(s, tab); |

258 | bb6f5690 | Fabrice Bellard | } |

259 | } |
||

260 | duration = gettime() - time_start; |
||

261 | if (duration >= 1000000) |
||

262 | ```
break;
``` |
||

263 | ```
nb_its *= 2;
``` |
||

264 | } |
||

265 | ```
printf("time: %0.1f us/transform [total time=%0.2f s its=%d]\n",
``` |
||

266 | ```
(double)duration / nb_its,
``` |
||

267 | (double)duration / 1000000.0, |
||

268 | nb_its); |
||

269 | } |
||

270 | |||

271 | ```
if (do_mdct) {
``` |
||

272 | ea0f841a | Fabrice Bellard | ff_mdct_end(m); |

273 | bb6f5690 | Fabrice Bellard | ```
} else {
``` |

274 | 68951ecf | Gildas Bazin | ff_fft_end(s); |

275 | bb6f5690 | Fabrice Bellard | } |

276 | return 0; |
||

277 | } |