Fix for an assertion when running the fixed point tests
[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 #include <stdio.h>
33
34 #include "mdct.h"
35 #include "stack_alloc.h"
36 #include "kiss_fft.h"
37 #include "mdct.h"
38 #include "modes.h"
39
40 #ifndef M_PI
41 #define M_PI 3.141592653
42 #endif
43
44 int ret = 0;
45 void check(kiss_fft_scalar  * in,kiss_fft_scalar  * out,int nfft,int isinverse)
46 {
47     int bin,k;
48     double errpow=0,sigpow=0;
49     double snr;
50     for (bin=0;bin<nfft/2;++bin) {
51         double ansr = 0;
52         double difr;
53
54         for (k=0;k<nfft;++k) {
55            double phase = 2*M_PI*(k+.5+.25*nfft)*(bin+.5)/nfft;
56            double re = cos(phase);
57
58            re /= nfft/4;
59
60            ansr += in[k] * re;
61         }
62         /*printf ("%f %f\n", ansr, out[bin]);*/
63         difr = ansr - out[bin];
64         errpow += difr*difr;
65         sigpow += ansr*ansr;
66     }
67     snr = 10*log10(sigpow/errpow);
68     printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
69     if (snr<60) {
70        printf( "** poor snr: %f **\n", snr);
71        ret = 1;
72     }
73 }
74
75 void check_inv(kiss_fft_scalar  * in,kiss_fft_scalar  * out,int nfft,int isinverse)
76 {
77    int bin,k;
78    double errpow=0,sigpow=0;
79    double snr;
80    for (bin=0;bin<nfft;++bin) {
81       double ansr = 0;
82       double difr;
83
84       for (k=0;k<nfft/2;++k) {
85          double phase = 2*M_PI*(bin+.5+.25*nfft)*(k+.5)/nfft;
86          double re = cos(phase);
87
88          /*re *= 2;*/
89
90          ansr += in[k] * re;
91       }
92       /*printf ("%f %f\n", ansr, out[bin]);*/
93       difr = ansr - out[bin];
94       errpow += difr*difr;
95       sigpow += ansr*ansr;
96    }
97    snr = 10*log10(sigpow/errpow);
98    printf("nfft=%d inverse=%d,snr = %f\n",nfft,isinverse,snr );
99    if (snr<60) {
100       printf( "** poor snr: %f **\n", snr);
101       ret = 1;
102    }
103 }
104
105
106 void test1d(int nfft,int isinverse,int arch)
107 {
108     size_t buflen = sizeof(kiss_fft_scalar)*nfft;
109     kiss_fft_scalar *in;
110     kiss_fft_scalar *in_copy;
111     kiss_fft_scalar *out;
112     opus_val16 *window;
113     int k;
114
115 #ifdef CUSTOM_MODES
116     int shift = 0;
117     const mdct_lookup *cfg;
118     mdct_lookup _cfg;
119     clt_mdct_init(&_cfg, nfft, 0, arch);
120     cfg = &_cfg;
121 #else
122     int shift;
123     const mdct_lookup *cfg;
124     CELTMode *mode = opus_custom_mode_create(48000, 960, NULL);
125     if (nfft == 1920) shift = 0;
126     else if (nfft == 960) shift = 1;
127     else if (nfft == 480) shift = 2;
128     else if (nfft == 240) shift = 3;
129     else return;
130     cfg = &mode->mdct;
131 #endif
132
133     in = (kiss_fft_scalar*)malloc(buflen);
134     in_copy = (kiss_fft_scalar*)malloc(buflen);
135     out = (kiss_fft_scalar*)malloc(buflen);
136     window = (opus_val16*)malloc(sizeof(opus_val16)*nfft/2);
137
138     for (k=0;k<nfft;++k) {
139         in[k] = (rand() % 32768) - 16384;
140     }
141
142     for (k=0;k<nfft/2;++k) {
143        window[k] = Q15ONE;
144     }
145     for (k=0;k<nfft;++k) {
146        in[k] *= 32768;
147     }
148
149     if (isinverse)
150     {
151        for (k=0;k<nfft;++k) {
152           in[k] /= nfft;
153        }
154     }
155
156     for (k=0;k<nfft;++k)
157        in_copy[k] = in[k];
158     /*for (k=0;k<nfft;++k) printf("%d %d ", in[k].r, in[k].i);printf("\n");*/
159
160     if (isinverse)
161     {
162        for (k=0;k<nfft;++k)
163           out[k] = 0;
164        clt_mdct_backward(cfg,in,out, window, nfft/2, shift, 1, arch);
165        /* apply TDAC because clt_mdct_backward() no longer does that */
166        for (k=0;k<nfft/4;++k)
167           out[nfft-k-1] = out[nfft/2+k];
168        check_inv(in,out,nfft,isinverse);
169     } else {
170        clt_mdct_forward(cfg,in,out,window, nfft/2, shift, 1, arch);
171        check(in_copy,out,nfft,isinverse);
172     }
173     /*for (k=0;k<nfft;++k) printf("%d %d ", out[k].r, out[k].i);printf("\n");*/
174
175
176     free(in);
177     free(in_copy);
178     free(out);
179     free(window);
180 #ifdef CUSTOM_MODES
181     clt_mdct_clear(&_cfg, arch);
182 #endif
183 }
184
185 int main(int argc,char ** argv)
186 {
187     ALLOC_STACK;
188     int arch = opus_select_arch();
189
190     if (argc>1) {
191         int k;
192         for (k=1;k<argc;++k) {
193             test1d(atoi(argv[k]),0,arch);
194             test1d(atoi(argv[k]),1,arch);
195         }
196     }else{
197         test1d(32,0,arch);
198         test1d(32,1,arch);
199         test1d(256,0,arch);
200         test1d(256,1,arch);
201         test1d(512,0,arch);
202         test1d(512,1,arch);
203         test1d(1024,0,arch);
204         test1d(1024,1,arch);
205         test1d(2048,0,arch);
206         test1d(2048,1,arch);
207 #ifndef RADIX_TWO_ONLY
208         test1d(36,0,arch);
209         test1d(36,1,arch);
210         test1d(40,0,arch);
211         test1d(40,1,arch);
212         test1d(60,0,arch);
213         test1d(60,1,arch);
214         test1d(120,0,arch);
215         test1d(120,1,arch);
216         test1d(240,0,arch);
217         test1d(240,1,arch);
218         test1d(480,0,arch);
219         test1d(480,1,arch);
220         test1d(960,0,arch);
221         test1d(960,1,arch);
222         test1d(1920,0,arch);
223         test1d(1920,1,arch);
224 #endif
225     }
226     return ret;
227 }