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