MDCT is in fixed-point now
[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(kiss_fft_scalar  * in,kiss_fft_scalar  * 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(kiss_fft_scalar  * in,kiss_fft_scalar  * 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(kiss_fft_scalar)*nfft;
74
75     kiss_fft_scalar  * in = (kiss_fft_scalar*)malloc(buflen);
76     kiss_fft_scalar  * out= (kiss_fft_scalar*)malloc(buflen);
77     int k;
78
79     mdct_init(&cfg, nfft);
80     for (k=0;k<nfft;++k) {
81         in[k] = (rand() % 32768) - 16384;
82     }
83
84 #ifdef DOUBLE_PRECISION
85     for (k=0;k<nfft;++k) {
86        in[k] *= 65536;
87     }
88 #endif
89     
90     if (isinverse)
91     {
92        for (k=0;k<nfft;++k) {
93           in[k] /= nfft;
94        }
95     }
96     
97     /*for (k=0;k<nfft;++k) printf("%d %d ", in[k].r, in[k].i);printf("\n");*/
98        
99     if (isinverse)
100     {
101        mdct_backward(&cfg,in,out);
102        check_inv(in,out,nfft,isinverse);
103     } else {
104        mdct_forward(&cfg,in,out);
105        check(in,out,nfft,isinverse);
106     }
107     /*for (k=0;k<nfft;++k) printf("%d %d ", out[k].r, out[k].i);printf("\n");*/
108
109
110     free(in);
111     free(out);
112     mdct_clear(&cfg);
113 }
114
115 int main(int argc,char ** argv)
116 {
117     if (argc>1) {
118         int k;
119         for (k=1;k<argc;++k) {
120             test1d(atoi(argv[k]),0);
121             test1d(atoi(argv[k]),1);
122         }
123     }else{
124         test1d(32,0);
125         test1d(32,1);
126         test1d(40,0);
127         test1d(40,1);
128         test1d(56,0);
129         test1d(56,1);
130         test1d(120,0);
131         test1d(120,1);
132         test1d(240,0);
133         test1d(240,1);
134         test1d(256,0);
135         test1d(256,1);
136         test1d(480,0);
137         test1d(480,1);
138         test1d(512,0);
139         test1d(512,1);
140     }
141     return ret;
142 }