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