Cleanup, de-inlining some math functions
[opus.git] / tests / dft-test.c
1 #ifdef HAVE_CONFIG_H
2 #include "config.h"
3 #endif
4
5 #include <stdio.h>
6 #include "kiss_fft.h"
7
8 #define CELT_C 
9 #include "../libcelt/stack_alloc.h"
10 #include "../libcelt/kiss_fft.c"
11 #include "../libcelt/mathops.c"
12 #include "../libcelt/entcode.c"
13
14
15 #ifndef M_PI
16 #define M_PI 3.141592653
17 #endif
18
19 #ifdef FIXED_DEBUG  
20 long long celt_mips=0;
21 #endif
22 int ret = 0;
23
24 void check(kiss_fft_cpx  * in,kiss_fft_cpx  * out,int nfft,int isinverse)
25 {
26     int bin,k;
27     double errpow=0,sigpow=0, snr;
28     
29     for (bin=0;bin<nfft;++bin) {
30         double ansr = 0;
31         double ansi = 0;
32         double difr;
33         double difi;
34
35         for (k=0;k<nfft;++k) {
36             double phase = -2*M_PI*bin*k/nfft;
37             double re = cos(phase);
38             double im = sin(phase);
39             if (isinverse)
40                 im = -im;
41
42             if (!isinverse)
43             {
44                re /= nfft;
45                im /= nfft;
46             }
47
48             ansr += in[k].r * re - in[k].i * im;
49             ansi += in[k].r * im + in[k].i * re;
50         }
51         /*printf ("%d %d ", (int)ansr, (int)ansi);*/
52         difr = ansr - out[bin].r;
53         difi = ansi - out[bin].i;
54         errpow += difr*difr + difi*difi;
55         sigpow += ansr*ansr+ansi*ansi;
56     }
57     snr = 10*log10(sigpow/errpow);
58     printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
59     if (snr<60) {
60        printf( "** poor snr: %f ** \n", snr);
61        ret = 1;
62     }
63 }
64
65 void test1d(int nfft,int isinverse)
66 {
67     size_t buflen = sizeof(kiss_fft_cpx)*nfft;
68
69     kiss_fft_cpx  * in = (kiss_fft_cpx*)malloc(buflen);
70     kiss_fft_cpx  * out= (kiss_fft_cpx*)malloc(buflen);
71     kiss_fft_cfg  cfg = kiss_fft_alloc(nfft,0,0);
72     int k;
73
74     for (k=0;k<nfft;++k) {
75         in[k].r = (rand() % 32767) - 16384;
76         in[k].i = (rand() % 32767) - 16384;
77     }
78
79 #ifdef DOUBLE_PRECISION
80     for (k=0;k<nfft;++k) {
81        in[k].r *= 32768;
82        in[k].i *= 32768;
83     }
84 #endif
85     
86     if (isinverse)
87     {
88        for (k=0;k<nfft;++k) {
89           in[k].r /= nfft;
90           in[k].i /= nfft;
91        }
92     }
93     
94     /*for (k=0;k<nfft;++k) printf("%d %d ", in[k].r, in[k].i);printf("\n");*/
95        
96     if (isinverse)
97        kiss_ifft(cfg,in,out);
98     else
99        kiss_fft(cfg,in,out);
100
101     /*for (k=0;k<nfft;++k) printf("%d %d ", out[k].r, out[k].i);printf("\n");*/
102
103     check(in,out,nfft,isinverse);
104
105     free(in);
106     free(out);
107     free(cfg);
108 }
109
110 int main(int argc,char ** argv)
111 {
112     ALLOC_STACK;
113     if (argc>1) {
114         int k;
115         for (k=1;k<argc;++k) {
116             test1d(atoi(argv[k]),0);
117             test1d(atoi(argv[k]),1);
118         }
119     }else{
120         test1d(32,0);
121         test1d(32,1);
122         test1d(128,0);
123         test1d(128,1);
124         test1d(256,0);
125         test1d(256,1);
126 #ifndef RADIX_TWO_ONLY
127         test1d(36,0);
128         test1d(36,1);
129         test1d(50,0);
130         test1d(50,1);
131         test1d(120,0);
132         test1d(120,1);
133 #endif
134     }
135     return ret;
136 }