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