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