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;
check(kiss_fft_scalar * in,kiss_fft_scalar * out,int nfft,int isinverse)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
check_inv(kiss_fft_scalar * in,kiss_fft_scalar * out,int nfft,int isinverse)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
test1d(int nfft,int isinverse)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
main(int argc,char ** argv)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 test1d(1024,0);
187 test1d(1024,1);
188 test1d(2048,0);
189 test1d(2048,1);
190 #ifndef RADIX_TWO_ONLY
191 test1d(36,0);
192 test1d(36,1);
193 test1d(40,0);
194 test1d(40,1);
195 test1d(60,0);
196 test1d(60,1);
197 test1d(120,0);
198 test1d(120,1);
199 test1d(240,0);
200 test1d(240,1);
201 test1d(480,0);
202 test1d(480,1);
203 test1d(960,0);
204 test1d(960,1);
205 test1d(1920,0);
206 test1d(1920,1);
207 #endif
208 }
209 return ret;
210 }
211