← Back to branch summary

~ubuntu-branches/ubuntu/hardy/avidemux/hardy

« back to all changes in this revision

Viewing changes to avidemux/ADM_lavcodec/fft-test.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T Chen
  • Date: 2006-12-15 17:13:20 UTC
  • mfrom: (1.1.6 upstream)
  • Revision ID: james.westby@ubuntu.com-20061215171320-w79pvpehxx2fr217
Tags: 1:2.3.0-0.0ubuntu1
* Merge from debian-multimedia.org, remaining Ubuntu change:
  - desktop file,
  - no support for ccache and make -j.
* Closes Ubuntu: #69614.

Show diffs side-by-side

added added

removed removed

Lines of Context:
avidemux/ADM_lavcodec/fft-test.c
 
1
/**
 
2
 * @file fft-test.c
 
3
 * FFT and MDCT tests.
 
4
 */
 
5
 
 
6
#include "dsputil.h"
 
7
#include <math.h>
 
8
#include <unistd.h>
 
9
#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
int64_t gettime(void)
 
113
{
 
114
    struct timeval tv;
 
115
    gettimeofday(&tv,NULL);
 
116
    return (int64_t)tv.tv_sec * 1000000 + tv.tv_usec;
 
117
}
 
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
            av_log(NULL, AV_LOG_ERROR, "ERROR %d: %f %f\n",
 
126
                   i, tab1[i], tab2[i]);
 
127
        }
 
128
    }
 
129
}
 
130
 
 
131
 
 
132
void help(void)
 
133
{
 
134
    av_log(NULL, AV_LOG_INFO,"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
            av_log(NULL, AV_LOG_INFO,"IMDCT");
 
193
        else
 
194
            av_log(NULL, AV_LOG_INFO,"MDCT");
 
195
        ff_mdct_init(m, fft_nbits, do_inverse);
 
196
    } else {
 
197
        if (do_inverse)
 
198
            av_log(NULL, AV_LOG_INFO,"IFFT");
 
199
        else
 
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);
 
203
    }
 
204
    av_log(NULL, AV_LOG_INFO," %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
    av_log(NULL, AV_LOG_INFO,"Checking...\n");
 
215
 
 
216
    if (do_mdct) {
 
217
        if (do_inverse) {
 
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);
 
221
        } else {
 
222
            mdct_ref((float *)tab_ref, (float *)tab1, fft_size);
 
223
 
 
224
            ff_mdct_calc(m, tab2, (float *)tab1, tabtmp);
 
225
 
 
226
            check_diff((float *)tab_ref, tab2, fft_size / 2);
 
227
        }
 
228
    } else {
 
229
        memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
 
230
        ff_fft_permute(s, tab);
 
231
        ff_fft_calc(s, tab);
 
232
 
 
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
        int64_t time_start, duration;
 
241
        int nb_its;
 
242
 
 
243
        av_log(NULL, AV_LOG_INFO,"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
                        ff_imdct_calc(m, (float *)tab, (float *)tab1, tabtmp);
 
252
                    } else {
 
253
                        ff_mdct_calc(m, (float *)tab, (float *)tab1, tabtmp);
 
254
                    }
 
255
                } else {
 
256
                    memcpy(tab, tab1, fft_size * sizeof(FFTComplex));
 
257
                    ff_fft_calc(s, tab);
 
258
                }
 
259
            }
 
260
            duration = gettime() - time_start;
 
261
            if (duration >= 1000000)
 
262
                break;
 
263
            nb_its *= 2;
 
264
        }
 
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,
 
268
               nb_its);
 
269
    }
 
270
 
 
271
    if (do_mdct) {
 
272
        ff_mdct_end(m);
 
273
    } else {
 
274
        ff_fft_end(s);
 
275
    }
 
276
    return 0;
 
277
}