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