MDCT now scales down by N/2 instead of N/4. The factor two is moved to the
[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 int ret = 0;
9 void check(float  * in,float  * out,int nfft,int isinverse)
10 {
11     int bin,k;
12     double errpow=0,sigpow=0;
13     double snr;
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/2;
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     snr = 10*log10(sigpow/errpow);
32     printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
33     if (snr<60) {
34        printf( "** poor snr: %f **\n", snr);
35        ret = 1;
36     }
37 }
38
39 void check_inv(float  * in,float  * out,int nfft,int isinverse)
40 {
41    int bin,k;
42    double errpow=0,sigpow=0;
43    double snr;
44    for (bin=0;bin<nfft;++bin) {
45       double ansr = 0;
46       double difr;
47
48       for (k=0;k<nfft/2;++k) {
49          double phase = 2*M_PI*(bin+.5+.25*nfft)*(k+.5)/nfft;
50          double re = cos(phase);
51
52          //re *= 2;
53
54          ansr += in[k] * re;
55       }
56       /*printf ("%f %f\n", ansr, out[bin]);*/
57       difr = ansr - out[bin];
58       errpow += difr*difr;
59       sigpow += ansr*ansr;
60    }
61    snr = 10*log10(sigpow/errpow);
62    printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
63    if (snr<60) {
64       printf( "** poor snr: %f **\n", snr);
65       ret = 1;
66    }
67 }
68
69
70 void test1d(int nfft,int isinverse)
71 {
72     mdct_lookup cfg;
73     size_t buflen = sizeof(float)*nfft;
74
75     float  * in = (float*)malloc(buflen);
76     float  * out= (float*)malloc(buflen);
77     int k;
78
79     mdct_init(&cfg, nfft);
80     for (k=0;k<nfft;++k) {
81         in[k] = (rand() % 65536) - 32768;
82     }
83
84 #ifdef DOUBLE_PRECISION
85     for (k=0;k<nfft;++k) {
86        in[k].r *= 65536;
87        in[k].i *= 65536;
88     }
89 #endif
90     
91     if (isinverse)
92     {
93        for (k=0;k<nfft;++k) {
94           in[k] /= nfft;
95        }
96     }
97     
98     /*for (k=0;k<nfft;++k) printf("%d %d ", in[k].r, in[k].i);printf("\n");*/
99        
100     if (isinverse)
101     {
102        mdct_backward(&cfg,in,out);
103        check_inv(in,out,nfft,isinverse);
104     } else {
105        mdct_forward(&cfg,in,out);
106        check(in,out,nfft,isinverse);
107     }
108     /*for (k=0;k<nfft;++k) printf("%d %d ", out[k].r, out[k].i);printf("\n");*/
109
110
111     free(in);
112     free(out);
113     mdct_clear(&cfg);
114 }
115
116 int main(int argc,char ** argv)
117 {
118     if (argc>1) {
119         int k;
120         for (k=1;k<argc;++k) {
121             test1d(atoi(argv[k]),0);
122             test1d(atoi(argv[k]),1);
123         }
124     }else{
125         test1d(32,0);
126         test1d(32,1);
127         test1d(40,0);
128         test1d(40,1);
129         test1d(56,0);
130         test1d(56,1);
131         test1d(120,0);
132         test1d(120,1);
133         test1d(240,0);
134         test1d(240,1);
135         test1d(256,0);
136         test1d(256,1);
137         test1d(480,0);
138         test1d(480,1);
139         test1d(512,0);
140         test1d(512,1);
141     }
142     return ret;
143 }