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 } 228
