Merge branch 'exp_api_change'
[opus.git] / tests / mdct-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 "mdct.h"
13 #define CELT_C 
14 #include "../libcelt/stack_alloc.h"
15
16 #include "../libcelt/kiss_fft.c"
17 #include "../libcelt/mdct.c"
18 #include "../libcelt/mathops.c"
19 #include "../libcelt/entcode.c"
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 void check(kiss_fft_scalar  * in,kiss_fft_scalar  * out,int nfft,int isinverse)
30 {
31     int bin,k;
32     double errpow=0,sigpow=0;
33     double snr;
34     for (bin=0;bin<nfft/2;++bin) {
35         double ansr = 0;
36         double difr;
37
38         for (k=0;k<nfft;++k) {
39            double phase = 2*M_PI*(k+.5+.25*nfft)*(bin+.5)/nfft;
40            double re = cos(phase);
41             
42            re /= nfft/4;
43
44            ansr += in[k] * re;
45         }
46         /*printf ("%f %f\n", ansr, out[bin]);*/
47         difr = ansr - out[bin];
48         errpow += difr*difr;
49         sigpow += ansr*ansr;
50     }
51     snr = 10*log10(sigpow/errpow);
52     printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
53     if (snr<60) {
54        printf( "** poor snr: %f **\n", snr);
55        ret = 1;
56     }
57 }
58
59 void check_inv(kiss_fft_scalar  * in,kiss_fft_scalar  * out,int nfft,int isinverse)
60 {
61    int bin,k;
62    double errpow=0,sigpow=0;
63    double snr;
64    for (bin=0;bin<nfft;++bin) {
65       double ansr = 0;
66       double difr;
67
68       for (k=0;k<nfft/2;++k) {
69          double phase = 2*M_PI*(bin+.5+.25*nfft)*(k+.5)/nfft;
70          double re = cos(phase);
71
72          //re *= 2;
73
74          ansr += in[k] * re;
75       }
76       /*printf ("%f %f\n", ansr, out[bin]);*/
77       difr = ansr - out[bin];
78       errpow += difr*difr;
79       sigpow += ansr*ansr;
80    }
81    snr = 10*log10(sigpow/errpow);
82    printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
83    if (snr<60) {
84       printf( "** poor snr: %f **\n", snr);
85       ret = 1;
86    }
87 }
88
89
90 void test1d(int nfft,int isinverse)
91 {
92     mdct_lookup cfg;
93     size_t buflen = sizeof(kiss_fft_scalar)*nfft;
94
95     kiss_fft_scalar  * in = (kiss_fft_scalar*)malloc(buflen);
96     kiss_fft_scalar  * out= (kiss_fft_scalar*)malloc(buflen);
97     celt_word16  * window= (celt_word16*)malloc(sizeof(celt_word16)*nfft/2);
98     int k;
99
100     clt_mdct_init(&cfg, nfft, 0);
101     for (k=0;k<nfft;++k) {
102         in[k] = (rand() % 32768) - 16384;
103     }
104
105     for (k=0;k<nfft/2;++k) {
106        window[k] = Q15ONE;
107     }
108 #ifdef DOUBLE_PRECISION
109     for (k=0;k<nfft;++k) {
110        in[k] *= 32768;
111     }
112 #endif
113     
114     if (isinverse)
115     {
116        for (k=0;k<nfft;++k) {
117           in[k] /= nfft;
118        }
119     }
120     
121     /*for (k=0;k<nfft;++k) printf("%d %d ", in[k].r, in[k].i);printf("\n");*/
122        
123     if (isinverse)
124     {
125        for (k=0;k<nfft;++k)
126           out[k] = 0;
127        clt_mdct_backward(&cfg,in,out, window, nfft/2, 0);
128        check_inv(in,out,nfft,isinverse);
129     } else {
130        clt_mdct_forward(&cfg,in,out,window, nfft/2, 0);
131        check(in,out,nfft,isinverse);
132     }
133     /*for (k=0;k<nfft;++k) printf("%d %d ", out[k].r, out[k].i);printf("\n");*/
134
135
136     free(in);
137     free(out);
138     clt_mdct_clear(&cfg);
139 }
140
141 int main(int argc,char ** argv)
142 {
143     ALLOC_STACK;
144     if (argc>1) {
145         int k;
146         for (k=1;k<argc;++k) {
147             test1d(atoi(argv[k]),0);
148             test1d(atoi(argv[k]),1);
149         }
150     }else{
151         test1d(32,0);
152         test1d(32,1);
153         test1d(256,0);
154         test1d(256,1);
155         test1d(512,0);
156         test1d(512,1);
157 #ifndef RADIX_TWO_ONLY
158         test1d(40,0);
159         test1d(40,1);
160         test1d(120,0);
161         test1d(120,1);
162         test1d(240,0);
163         test1d(240,1);
164         test1d(480,0);
165         test1d(480,1);
166 #endif
167     }
168     return ret;
169 }