testcase for the MDCT and IMDCT
[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
8
9 void check(float  * in,float  * out,int nfft,int isinverse)
10 {
11     int bin,k;
12     double errpow=0,sigpow=0;
13     
14     for (bin=0;bin<nfft/2;++bin) {
15         double ansr = 0;
16         double difr;
17
18         for (k=0;k<nfft;++k) {
19            double phase = 2*M_PI*(k+.5+.25*nfft)*(bin+.5)/nfft;
20            double re = cos(phase);
21             
22            re /= nfft/4;
23
24            ansr += in[k] * re;
25         }
26         /*printf ("%f %f\n", ansr, out[bin]);*/
27         difr = ansr - out[bin];
28         errpow += difr*difr;
29         sigpow += ansr*ansr;
30     }
31     printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,10*log10(sigpow/errpow) );
32 }
33
34 void check_inv(float  * in,float  * out,int nfft,int isinverse)
35 {
36    int bin,k;
37    double errpow=0,sigpow=0;
38     
39    for (bin=0;bin<nfft;++bin) {
40       double ansr = 0;
41       double difr;
42
43       for (k=0;k<nfft/2;++k) {
44          double phase = 2*M_PI*(bin+.5+.25*nfft)*(k+.5)/nfft;
45          double re = cos(phase);
46             
47          ansr += in[k] * re;
48       }
49       /*printf ("%f %f\n", ansr, out[bin]);*/
50       difr = ansr - out[bin];
51       errpow += difr*difr;
52       sigpow += ansr*ansr;
53    }
54    printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,10*log10(sigpow/errpow) );
55 }
56
57
58 void test1d(int nfft,int isinverse)
59 {
60     mdct_lookup cfg;
61     size_t buflen = sizeof(float)*nfft;
62
63     float  * in = (float*)malloc(buflen);
64     float  * out= (float*)malloc(buflen);
65     int k;
66
67     mdct_init(&cfg, nfft);
68     for (k=0;k<nfft;++k) {
69         in[k] = (rand() % 65536) - 32768;
70     }
71
72 #ifdef DOUBLE_PRECISION
73     for (k=0;k<nfft;++k) {
74        in[k].r *= 65536;
75        in[k].i *= 65536;
76     }
77 #endif
78     
79     if (isinverse)
80     {
81        for (k=0;k<nfft;++k) {
82           in[k] /= nfft;
83        }
84     }
85     
86     /*for (k=0;k<nfft;++k) printf("%d %d ", in[k].r, in[k].i);printf("\n");*/
87        
88     if (isinverse)
89     {
90        mdct_backward(&cfg,in,out);
91        check_inv(in,out,nfft,isinverse);
92     } else {
93        mdct_forward(&cfg,in,out);
94        check(in,out,nfft,isinverse);
95     }
96     /*for (k=0;k<nfft;++k) printf("%d %d ", out[k].r, out[k].i);printf("\n");*/
97
98
99     free(in);
100     free(out);
101     mdct_clear(&cfg);
102 }
103
104 int main(int argc,char ** argv)
105 {
106     if (argc>1) {
107         int k;
108         for (k=1;k<argc;++k) {
109             test1d(atoi(argv[k]),0);
110             test1d(atoi(argv[k]),1);
111         }
112     }else{
113         test1d(32,0);
114         test1d(32,1);
115         exit(0);
116         test1d(36,0);
117         test1d(36,1);
118         test1d(50,0);
119         test1d(50,1);
120         test1d(120,0);
121         test1d(120,1);
122     }
123     return 0;
124 }