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;
check(kiss_fft_scalar * in,kiss_fft_scalar * out,int nfft,int isinverse)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
check_inv(kiss_fft_scalar * in,kiss_fft_scalar * out,int nfft,int isinverse)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
test1d(int nfft,int isinverse,int arch)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
main(int argc,char ** argv)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 }
228