1 /* Copyright (c) 2011-2012 Xiph.Org Foundation, Mozilla Corporation 2 Written by Jean-Marc Valin and Timothy B. Terriberry */ 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 #include <stdio.h> 29 #include <stdlib.h> 30 #include <math.h> 31 #include <string.h> 32 33 #define OPUS_PI (3.14159265F) 34 35 #define OPUS_COSF(_x) ((float)cos(_x)) 36 #define OPUS_SINF(_x) ((float)sin(_x)) 37 38 static void *check_alloc(void *_ptr){ 39 if(_ptr==NULL){ 40 fprintf(stderr,"Out of memory.\n"); 41 exit(EXIT_FAILURE); 42 } 43 return _ptr; 44 } 45 46 static void *opus_malloc(size_t _size){ 47 return check_alloc(malloc(_size)); 48 } 49 50 static void *opus_realloc(void *_ptr,size_t _size){ 51 return check_alloc(realloc(_ptr,_size)); 52 } 53 54 static size_t read_pcm16(float **_samples,FILE *_fin,int _nchannels){ 55 unsigned char buf[1024]; 56 float *samples; 57 size_t nsamples; 58 size_t csamples; 59 size_t xi; 60 size_t nread; 61 samples=NULL; 62 nsamples=csamples=0; 63 for(;;){ 64 nread=fread(buf,2*_nchannels,1024/(2*_nchannels),_fin); 65 if(nread<=0)break; 66 if(nsamples+nread>csamples){ 67 do csamples=csamples<<1|1; 68 while(nsamples+nread>csamples); 69 samples=(float *)opus_realloc(samples, 70 _nchannels*csamples*sizeof(*samples)); 71 } 72 for(xi=0;xi<nread;xi++){ 73 int ci; 74 for(ci=0;ci<_nchannels;ci++){ 75 int s; 76 s=buf[2*(xi*_nchannels+ci)+1]<<8|buf[2*(xi*_nchannels+ci)]; 77 s=((s&0xFFFF)^0x8000)-0x8000; 78 samples[(nsamples+xi)*_nchannels+ci]=s; 79 } 80 } 81 nsamples+=nread; 82 } 83 *_samples=(float *)opus_realloc(samples, 84 _nchannels*nsamples*sizeof(*samples)); 85 return nsamples; 86 } 87 88 static void band_energy(float *_out,float *_ps,const int *_bands,int _nbands, 89 const float *_in,int _nchannels,size_t _nframes,int _window_sz, 90 int _step,int _downsample){ 91 float *window; 92 float *x; 93 float *c; 94 float *s; 95 size_t xi; 96 int xj; 97 int ps_sz; 98 window=(float *)opus_malloc((3+_nchannels)*_window_sz*sizeof(*window)); 99 c=window+_window_sz; 100 s=c+_window_sz; 101 x=s+_window_sz; 102 ps_sz=_window_sz/2; 103 for(xj=0;xj<_window_sz;xj++){ 104 window[xj]=0.5F-0.5F*OPUS_COSF((2*OPUS_PI/(_window_sz-1))*xj); 105 } 106 for(xj=0;xj<_window_sz;xj++){ 107 c[xj]=OPUS_COSF((2*OPUS_PI/_window_sz)*xj); 108 } 109 for(xj=0;xj<_window_sz;xj++){ 110 s[xj]=OPUS_SINF((2*OPUS_PI/_window_sz)*xj); 111 } 112 for(xi=0;xi<_nframes;xi++){ 113 int ci; 114 int xk; 115 int bi; 116 for(ci=0;ci<_nchannels;ci++){ 117 for(xk=0;xk<_window_sz;xk++){ 118 x[ci*_window_sz+xk]=window[xk]*_in[(xi*_step+xk)*_nchannels+ci]; 119 } 120 } 121 for(bi=xj=0;bi<_nbands;bi++){ 122 float p[2]={0}; 123 for(;xj<_bands[bi+1];xj++){ 124 for(ci=0;ci<_nchannels;ci++){ 125 float re; 126 float im; 127 int ti; 128 ti=0; 129 re=im=0; 130 for(xk=0;xk<_window_sz;xk++){ 131 re+=c[ti]*x[ci*_window_sz+xk]; 132 im-=s[ti]*x[ci*_window_sz+xk]; 133 ti+=xj; 134 if(ti>=_window_sz)ti-=_window_sz; 135 } 136 re*=_downsample; 137 im*=_downsample; 138 _ps[(xi*ps_sz+xj)*_nchannels+ci]=re*re+im*im+100000; 139 p[ci]+=_ps[(xi*ps_sz+xj)*_nchannels+ci]; 140 } 141 } 142 if(_out){ 143 _out[(xi*_nbands+bi)*_nchannels]=p[0]/(_bands[bi+1]-_bands[bi]); 144 if(_nchannels==2){ 145 _out[(xi*_nbands+bi)*_nchannels+1]=p[1]/(_bands[bi+1]-_bands[bi]); 146 } 147 } 148 } 149 } 150 free(window); 151 } 152 153 #define NBANDS (21) 154 #define NFREQS (240) 155 156 /*Bands on which we compute the pseudo-NMR (Bark-derived 157 CELT bands).*/ 158 static const int BANDS[NBANDS+1]={ 159 0,2,4,6,8,10,12,14,16,20,24,28,32,40,48,56,68,80,96,120,156,200 160 }; 161 162 #define TEST_WIN_SIZE (480) 163 #define TEST_WIN_STEP (120) 164 165 int main(int _argc,const char **_argv){ 166 FILE *fin1; 167 FILE *fin2; 168 float *x; 169 float *y; 170 float *xb; 171 float *X; 172 float *Y; 173 double err; 174 float Q; 175 size_t xlength; 176 size_t ylength; 177 size_t nframes; 178 size_t xi; 179 int ci; 180 int xj; 181 int bi; 182 int nchannels; 183 unsigned rate; 184 int downsample; 185 int ybands; 186 int yfreqs; 187 int max_compare; 188 if(_argc<3||_argc>6){ 189 fprintf(stderr,"Usage: %s [-s] [-r rate2] <file1.sw> <file2.sw>\n", 190 _argv[0]); 191 return EXIT_FAILURE; 192 } 193 nchannels=1; 194 if(strcmp(_argv[1],"-s")==0){ 195 nchannels=2; 196 _argv++; 197 } 198 rate=48000; 199 ybands=NBANDS; 200 yfreqs=NFREQS; 201 downsample=1; 202 if(strcmp(_argv[1],"-r")==0){ 203 rate=atoi(_argv[2]); 204 if(rate!=8000&&rate!=12000&&rate!=16000&&rate!=24000&&rate!=48000){ 205 fprintf(stderr, 206 "Sampling rate must be 8000, 12000, 16000, 24000, or 48000\n"); 207 return EXIT_FAILURE; 208 } 209 downsample=48000/rate; 210 switch(rate){ 211 case 8000:ybands=13;break; 212 case 12000:ybands=15;break; 213 case 16000:ybands=17;break; 214 case 24000:ybands=19;break; 215 } 216 yfreqs=NFREQS/downsample; 217 _argv+=2; 218 } 219 fin1=fopen(_argv[1],"rb"); 220 if(fin1==NULL){ 221 fprintf(stderr,"Error opening '%s'.\n",_argv[1]); 222 return EXIT_FAILURE; 223 } 224 fin2=fopen(_argv[2],"rb"); 225 if(fin2==NULL){ 226 fprintf(stderr,"Error opening '%s'.\n",_argv[2]); 227 fclose(fin1); 228 return EXIT_FAILURE; 229 } 230 /*Read in the data and allocate scratch space.*/ 231 xlength=read_pcm16(&x,fin1,2); 232 if(nchannels==1){ 233 for(xi=0;xi<xlength;xi++)x[xi]=.5*(x[2*xi]+x[2*xi+1]); 234 } 235 fclose(fin1); 236 ylength=read_pcm16(&y,fin2,nchannels); 237 fclose(fin2); 238 if(xlength!=ylength*downsample){ 239 fprintf(stderr,"Sample counts do not match (%lu!=%lu).\n", 240 (unsigned long)xlength,(unsigned long)ylength*downsample); 241 return EXIT_FAILURE; 242 } 243 if(xlength<TEST_WIN_SIZE){ 244 fprintf(stderr,"Insufficient sample data (%lu<%i).\n", 245 (unsigned long)xlength,TEST_WIN_SIZE); 246 return EXIT_FAILURE; 247 } 248 nframes=(xlength-TEST_WIN_SIZE+TEST_WIN_STEP)/TEST_WIN_STEP; 249 xb=(float *)opus_malloc(nframes*NBANDS*nchannels*sizeof(*xb)); 250 X=(float *)opus_malloc(nframes*NFREQS*nchannels*sizeof(*X)); 251 Y=(float *)opus_malloc(nframes*yfreqs*nchannels*sizeof(*Y)); 252 /*Compute the per-band spectral energy of the original signal 253 and the error.*/ 254 band_energy(xb,X,BANDS,NBANDS,x,nchannels,nframes, 255 TEST_WIN_SIZE,TEST_WIN_STEP,1); 256 free(x); 257 band_energy(NULL,Y,BANDS,ybands,y,nchannels,nframes, 258 TEST_WIN_SIZE/downsample,TEST_WIN_STEP/downsample,downsample); 259 free(y); 260 for(xi=0;xi<nframes;xi++){ 261 /*Frequency masking (low to high): 10 dB/Bark slope.*/ 262 for(bi=1;bi<NBANDS;bi++){ 263 for(ci=0;ci<nchannels;ci++){ 264 xb[(xi*NBANDS+bi)*nchannels+ci]+= 265 0.1F*xb[(xi*NBANDS+bi-1)*nchannels+ci]; 266 } 267 } 268 /*Frequency masking (high to low): 15 dB/Bark slope.*/ 269 for(bi=NBANDS-1;bi-->0;){ 270 for(ci=0;ci<nchannels;ci++){ 271 xb[(xi*NBANDS+bi)*nchannels+ci]+= 272 0.03F*xb[(xi*NBANDS+bi+1)*nchannels+ci]; 273 } 274 } 275 if(xi>0){ 276 /*Temporal masking: -3 dB/2.5ms slope.*/ 277 for(bi=0;bi<NBANDS;bi++){ 278 for(ci=0;ci<nchannels;ci++){ 279 xb[(xi*NBANDS+bi)*nchannels+ci]+= 280 0.5F*xb[((xi-1)*NBANDS+bi)*nchannels+ci]; 281 } 282 } 283 } 284 /* Allowing some cross-talk */ 285 if(nchannels==2){ 286 for(bi=0;bi<NBANDS;bi++){ 287 float l,r; 288 l=xb[(xi*NBANDS+bi)*nchannels+0]; 289 r=xb[(xi*NBANDS+bi)*nchannels+1]; 290 xb[(xi*NBANDS+bi)*nchannels+0]+=0.01F*r; 291 xb[(xi*NBANDS+bi)*nchannels+1]+=0.01F*l; 292 } 293 } 294 295 /* Apply masking */ 296 for(bi=0;bi<ybands;bi++){ 297 for(xj=BANDS[bi];xj<BANDS[bi+1];xj++){ 298 for(ci=0;ci<nchannels;ci++){ 299 X[(xi*NFREQS+xj)*nchannels+ci]+= 300 0.1F*xb[(xi*NBANDS+bi)*nchannels+ci]; 301 Y[(xi*yfreqs+xj)*nchannels+ci]+= 302 0.1F*xb[(xi*NBANDS+bi)*nchannels+ci]; 303 } 304 } 305 } 306 } 307 308 /* Average of consecutive frames to make comparison slightly less sensitive */ 309 for(bi=0;bi<ybands;bi++){ 310 for(xj=BANDS[bi];xj<BANDS[bi+1];xj++){ 311 for(ci=0;ci<nchannels;ci++){ 312 float xtmp; 313 float ytmp; 314 xtmp = X[xj*nchannels+ci]; 315 ytmp = Y[xj*nchannels+ci]; 316 for(xi=1;xi<nframes;xi++){ 317 float xtmp2; 318 float ytmp2; 319 xtmp2 = X[(xi*NFREQS+xj)*nchannels+ci]; 320 ytmp2 = Y[(xi*yfreqs+xj)*nchannels+ci]; 321 X[(xi*NFREQS+xj)*nchannels+ci] += xtmp; 322 Y[(xi*yfreqs+xj)*nchannels+ci] += ytmp; 323 xtmp = xtmp2; 324 ytmp = ytmp2; 325 } 326 } 327 } 328 } 329 330 /*If working at a lower sampling rate, don't take into account the last 331 300 Hz to allow for different transition bands. 332 For 12 kHz, we don't skip anything, because the last band already skips 333 400 Hz.*/ 334 if(rate==48000)max_compare=BANDS[NBANDS]; 335 else if(rate==12000)max_compare=BANDS[ybands]; 336 else max_compare=BANDS[ybands]-3; 337 err=0; 338 for(xi=0;xi<nframes;xi++){ 339 double Ef; 340 Ef=0; 341 for(bi=0;bi<ybands;bi++){ 342 double Eb; 343 Eb=0; 344 for(xj=BANDS[bi];xj<BANDS[bi+1]&&xj<max_compare;xj++){ 345 for(ci=0;ci<nchannels;ci++){ 346 float re; 347 float im; 348 re=Y[(xi*yfreqs+xj)*nchannels+ci]/X[(xi*NFREQS+xj)*nchannels+ci]; 349 im=re-log(re)-1; 350 /*Make comparison less sensitive around the SILK/CELT cross-over to 351 allow for mode freedom in the filters.*/ 352 if(xj>=79&&xj<=81)im*=0.1F; 353 if(xj==80)im*=0.1F; 354 Eb+=im; 355 } 356 } 357 Eb /= (BANDS[bi+1]-BANDS[bi])*nchannels; 358 Ef += Eb*Eb; 359 } 360 /*Using a fixed normalization value means we're willing to accept slightly 361 lower quality for lower sampling rates.*/ 362 Ef/=NBANDS; 363 Ef*=Ef; 364 err+=Ef*Ef; 365 } 366 free(xb); 367 free(X); 368 free(Y); 369 err=pow(err/nframes,1.0/16); 370 Q=100*(1-0.5*log(1+err)/log(1.13)); 371 if(Q<0){ 372 fprintf(stderr,"Test vector FAILS\n"); 373 fprintf(stderr,"Internal weighted error is %f\n",err); 374 return EXIT_FAILURE; 375 } 376 else{ 377 fprintf(stderr,"Test vector PASSES\n"); 378 fprintf(stderr, 379 "Opus quality metric: %.1f %% (internal weighted error is %f)\n",Q,err); 380 return EXIT_SUCCESS; 381 } 382 } 383
