Updates MDCT unit test
[opus.git] / celt / tests / test_unit_mdct.c
1 /* Copyright (c) 2008-2011 Xiph.Org Foundation
2    Written by Jean-Marc Valin */
3 /*
4    Redistribution and use in source and binary forms, with or without
5    modification, are permitted provided that the following conditions
6    are met:
7
8    - Redistributions of source code must retain the above copyright
9    notice, this list of conditions and the following disclaimer.
10
11    - Redistributions in binary form must reproduce the above copyright
12    notice, this list of conditions and the following disclaimer in the
13    documentation and/or other materials provided with the distribution.
14
15    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
19    OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
20    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
21    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
22    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
23    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
24    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 */
27
28 #ifdef HAVE_CONFIG_H
29 #include "config.h"
30 #endif
31
32 #define SKIP_CONFIG_H
33
34 #ifndef CUSTOM_MODES
35 #define CUSTOM_MODES
36 #endif
37
38 #include <stdio.h>
39
40 #define CELT_C
41 #include "mdct.h"
42 #include "stack_alloc.h"
43
44 #include "kiss_fft.c"
45 #include "mdct.c"
46 #include "mathops.c"
47 #include "entcode.c"
48
49 #ifndef M_PI
50 #define M_PI 3.141592653
51 #endif
52
53 int ret = 0;
54 void check(kiss_fft_scalar  * in,kiss_fft_scalar  * out,int nfft,int isinverse)
55 {
56     int bin,k;
57     double errpow=0,sigpow=0;
58     double snr;
59     for (bin=0;bin<nfft/2;++bin) {
60         double ansr = 0;
61         double difr;
62
63         for (k=0;k<nfft;++k) {
64            double phase = 2*M_PI*(k+.5+.25*nfft)*(bin+.5)/nfft;
65            double re = cos(phase);
66
67            re /= nfft/4;
68
69            ansr += in[k] * re;
70         }
71         /*printf ("%f %f\n", ansr, out[bin]);*/
72         difr = ansr - out[bin];
73         errpow += difr*difr;
74         sigpow += ansr*ansr;
75     }
76     snr = 10*log10(sigpow/errpow);
77     printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
78     if (snr<60) {
79        printf( "** poor snr: %f **\n", snr);
80        ret = 1;
81     }
82 }
83
84 void check_inv(kiss_fft_scalar  * in,kiss_fft_scalar  * out,int nfft,int isinverse)
85 {
86    int bin,k;
87    double errpow=0,sigpow=0;
88    double snr;
89    for (bin=0;bin<nfft;++bin) {
90       double ansr = 0;
91       double difr;
92
93       for (k=0;k<nfft/2;++k) {
94          double phase = 2*M_PI*(bin+.5+.25*nfft)*(k+.5)/nfft;
95          double re = cos(phase);
96
97          /*re *= 2;*/
98
99          ansr += in[k] * re;
100       }
101       /*printf ("%f %f\n", ansr, out[bin]);*/
102       difr = ansr - out[bin];
103       errpow += difr*difr;
104       sigpow += ansr*ansr;
105    }
106    snr = 10*log10(sigpow/errpow);
107    printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
108    if (snr<60) {
109       printf( "** poor snr: %f **\n", snr);
110       ret = 1;
111    }
112 }
113
114
115 void test1d(int nfft,int isinverse)
116 {
117     mdct_lookup cfg;
118     size_t buflen = sizeof(kiss_fft_scalar)*nfft;
119
120     kiss_fft_scalar  * in = (kiss_fft_scalar*)malloc(buflen);
121     kiss_fft_scalar  * in_copy = (kiss_fft_scalar*)malloc(buflen);
122     kiss_fft_scalar  * out= (kiss_fft_scalar*)malloc(buflen);
123     opus_val16  * window= (opus_val16*)malloc(sizeof(opus_val16)*nfft/2);
124     int k;
125
126     clt_mdct_init(&cfg, nfft, 0);
127     for (k=0;k<nfft;++k) {
128         in[k] = (rand() % 32768) - 16384;
129     }
130
131     for (k=0;k<nfft/2;++k) {
132        window[k] = Q15ONE;
133     }
134     for (k=0;k<nfft;++k) {
135        in[k] *= 32768;
136     }
137
138     if (isinverse)
139     {
140        for (k=0;k<nfft;++k) {
141           in[k] /= nfft;
142        }
143     }
144
145     for (k=0;k<nfft;++k)
146        in_copy[k] = in[k];
147     /*for (k=0;k<nfft;++k) printf("%d %d ", in[k].r, in[k].i);printf("\n");*/
148
149     if (isinverse)
150     {
151        for (k=0;k<nfft;++k)
152           out[k] = 0;
153        clt_mdct_backward(&cfg,in,out, window, nfft/2, 0, 1);
154        /* apply TDAC because clt_mdct_backward() no longer does that */
155        for (k=0;k<nfft/4;++k)
156           out[nfft-k-1] = out[nfft/2+k];
157        check_inv(in,out,nfft,isinverse);
158     } else {
159        clt_mdct_forward(&cfg,in,out,window, nfft/2, 0, 1);
160        check(in_copy,out,nfft,isinverse);
161     }
162     /*for (k=0;k<nfft;++k) printf("%d %d ", out[k].r, out[k].i);printf("\n");*/
163
164
165     free(in);
166     free(out);
167     clt_mdct_clear(&cfg);
168 }
169
170 int main(int argc,char ** argv)
171 {
172     ALLOC_STACK;
173     if (argc>1) {
174         int k;
175         for (k=1;k<argc;++k) {
176             test1d(atoi(argv[k]),0);
177             test1d(atoi(argv[k]),1);
178         }
179     }else{
180         test1d(32,0);
181         test1d(32,1);
182         test1d(256,0);
183         test1d(256,1);
184         test1d(512,0);
185         test1d(512,1);
186 #ifndef RADIX_TWO_ONLY
187         test1d(40,0);
188         test1d(40,1);
189         test1d(120,0);
190         test1d(120,1);
191         test1d(240,0);
192         test1d(240,1);
193         test1d(480,0);
194         test1d(480,1);
195 #endif
196     }
197     return ret;
198 }