Simplifying the stereo energy quantisation functions now that energy is
[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 int ret = 0;
18 void check(kiss_fft_scalar  * in,kiss_fft_scalar  * out,int nfft,int isinverse)
19 {
20     int bin,k;
21     double errpow=0,sigpow=0;
22     double snr;
23     for (bin=0;bin<nfft/2;++bin) {
24         double ansr = 0;
25         double difr;
26
27         for (k=0;k<nfft;++k) {
28            double phase = 2*M_PI*(k+.5+.25*nfft)*(bin+.5)/nfft;
29            double re = cos(phase);
30             
31            re /= nfft/4;
32
33            ansr += in[k] * re;
34         }
35         /*printf ("%f %f\n", ansr, out[bin]);*/
36         difr = ansr - out[bin];
37         errpow += difr*difr;
38         sigpow += ansr*ansr;
39     }
40     snr = 10*log10(sigpow/errpow);
41     printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
42     if (snr<60) {
43        printf( "** poor snr: %f **\n", snr);
44        ret = 1;
45     }
46 }
47
48 void check_inv(kiss_fft_scalar  * in,kiss_fft_scalar  * out,int nfft,int isinverse)
49 {
50    int bin,k;
51    double errpow=0,sigpow=0;
52    double snr;
53    for (bin=0;bin<nfft;++bin) {
54       double ansr = 0;
55       double difr;
56
57       for (k=0;k<nfft/2;++k) {
58          double phase = 2*M_PI*(bin+.5+.25*nfft)*(k+.5)/nfft;
59          double re = cos(phase);
60
61          //re *= 2;
62
63          ansr += in[k] * re;
64       }
65       /*printf ("%f %f\n", ansr, out[bin]);*/
66       difr = ansr - out[bin];
67       errpow += difr*difr;
68       sigpow += ansr*ansr;
69    }
70    snr = 10*log10(sigpow/errpow);
71    printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
72    if (snr<60) {
73       printf( "** poor snr: %f **\n", snr);
74       ret = 1;
75    }
76 }
77
78
79 void test1d(int nfft,int isinverse)
80 {
81     mdct_lookup cfg;
82     size_t buflen = sizeof(kiss_fft_scalar)*nfft;
83
84     kiss_fft_scalar  * in = (kiss_fft_scalar*)malloc(buflen);
85     kiss_fft_scalar  * out= (kiss_fft_scalar*)malloc(buflen);
86     celt_word16_t  * window= (celt_word16_t*)malloc(sizeof(celt_word16_t)*nfft/2);
87     int k;
88
89     mdct_init(&cfg, nfft);
90     for (k=0;k<nfft;++k) {
91         in[k] = (rand() % 32768) - 16384;
92     }
93
94     for (k=0;k<nfft/2;++k) {
95        window[k] = Q15ONE;
96     }
97 #ifdef DOUBLE_PRECISION
98     for (k=0;k<nfft;++k) {
99        in[k] *= 32768;
100     }
101 #endif
102     
103     if (isinverse)
104     {
105        for (k=0;k<nfft;++k) {
106           in[k] /= nfft;
107        }
108     }
109     
110     /*for (k=0;k<nfft;++k) printf("%d %d ", in[k].r, in[k].i);printf("\n");*/
111        
112     if (isinverse)
113     {
114        for (k=0;k<nfft;++k)
115           out[k] = 0;
116        mdct_backward(&cfg,in,out, window, nfft/2);
117        check_inv(in,out,nfft,isinverse);
118     } else {
119        mdct_forward(&cfg,in,out,window, nfft/2);
120        check(in,out,nfft,isinverse);
121     }
122     /*for (k=0;k<nfft;++k) printf("%d %d ", out[k].r, out[k].i);printf("\n");*/
123
124
125     free(in);
126     free(out);
127     mdct_clear(&cfg);
128 }
129
130 int main(int argc,char ** argv)
131 {
132     ALLOC_STACK;
133     if (argc>1) {
134         int k;
135         for (k=1;k<argc;++k) {
136             test1d(atoi(argv[k]),0);
137             test1d(atoi(argv[k]),1);
138         }
139     }else{
140         test1d(32,0);
141         test1d(32,1);
142         test1d(256,0);
143         test1d(256,1);
144         test1d(512,0);
145         test1d(512,1);
146 #ifndef RADIX_TWO_ONLY
147         test1d(40,0);
148         test1d(40,1);
149         test1d(56,0);
150         test1d(56,1);
151         test1d(120,0);
152         test1d(120,1);
153         test1d(240,0);
154         test1d(240,1);
155         test1d(480,0);
156         test1d(480,1);
157 #endif
158     }
159     return ret;
160 }