Alexandre Lision | 744f742 | 2013-09-25 11:39:37 -0400 | [diff] [blame] | 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 | err=pow(err/nframes,1.0/16); |
| 367 | Q=100*(1-0.5*log(1+err)/log(1.13)); |
| 368 | if(Q<0){ |
| 369 | fprintf(stderr,"Test vector FAILS\n"); |
| 370 | fprintf(stderr,"Internal weighted error is %f\n",err); |
| 371 | return EXIT_FAILURE; |
| 372 | } |
| 373 | else{ |
| 374 | fprintf(stderr,"Test vector PASSES\n"); |
| 375 | fprintf(stderr, |
| 376 | "Opus quality metric: %.1f %% (internal weighted error is %f)\n",Q,err); |
| 377 | return EXIT_SUCCESS; |
| 378 | } |
| 379 | } |