| /* Copyright (c) 2011-2012 Xiph.Org Foundation, Mozilla Corporation |
| Written by Jean-Marc Valin and Timothy B. Terriberry */ |
| /* |
| Redistribution and use in source and binary forms, with or without |
| modification, are permitted provided that the following conditions |
| are met: |
| |
| - Redistributions of source code must retain the above copyright |
| notice, this list of conditions and the following disclaimer. |
| |
| - Redistributions in binary form must reproduce the above copyright |
| notice, this list of conditions and the following disclaimer in the |
| documentation and/or other materials provided with the distribution. |
| |
| THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
| ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
| LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
| A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER |
| OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
| EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
| PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR |
| PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF |
| LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING |
| NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
| SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
| */ |
| |
| #include <stdio.h> |
| #include <stdlib.h> |
| #include <math.h> |
| #include <string.h> |
| |
| #define OPUS_PI (3.14159265F) |
| |
| #define OPUS_COSF(_x) ((float)cos(_x)) |
| #define OPUS_SINF(_x) ((float)sin(_x)) |
| |
| static void *check_alloc(void *_ptr){ |
| if(_ptr==NULL){ |
| fprintf(stderr,"Out of memory.\n"); |
| exit(EXIT_FAILURE); |
| } |
| return _ptr; |
| } |
| |
| static void *opus_malloc(size_t _size){ |
| return check_alloc(malloc(_size)); |
| } |
| |
| static void *opus_realloc(void *_ptr,size_t _size){ |
| return check_alloc(realloc(_ptr,_size)); |
| } |
| |
| static size_t read_pcm16(float **_samples,FILE *_fin,int _nchannels){ |
| unsigned char buf[1024]; |
| float *samples; |
| size_t nsamples; |
| size_t csamples; |
| size_t xi; |
| size_t nread; |
| samples=NULL; |
| nsamples=csamples=0; |
| for(;;){ |
| nread=fread(buf,2*_nchannels,1024/(2*_nchannels),_fin); |
| if(nread<=0)break; |
| if(nsamples+nread>csamples){ |
| do csamples=csamples<<1|1; |
| while(nsamples+nread>csamples); |
| samples=(float *)opus_realloc(samples, |
| _nchannels*csamples*sizeof(*samples)); |
| } |
| for(xi=0;xi<nread;xi++){ |
| int ci; |
| for(ci=0;ci<_nchannels;ci++){ |
| int s; |
| s=buf[2*(xi*_nchannels+ci)+1]<<8|buf[2*(xi*_nchannels+ci)]; |
| s=((s&0xFFFF)^0x8000)-0x8000; |
| samples[(nsamples+xi)*_nchannels+ci]=s; |
| } |
| } |
| nsamples+=nread; |
| } |
| *_samples=(float *)opus_realloc(samples, |
| _nchannels*nsamples*sizeof(*samples)); |
| return nsamples; |
| } |
| |
| static void band_energy(float *_out,float *_ps,const int *_bands,int _nbands, |
| const float *_in,int _nchannels,size_t _nframes,int _window_sz, |
| int _step,int _downsample){ |
| float *window; |
| float *x; |
| float *c; |
| float *s; |
| size_t xi; |
| int xj; |
| int ps_sz; |
| window=(float *)opus_malloc((3+_nchannels)*_window_sz*sizeof(*window)); |
| c=window+_window_sz; |
| s=c+_window_sz; |
| x=s+_window_sz; |
| ps_sz=_window_sz/2; |
| for(xj=0;xj<_window_sz;xj++){ |
| window[xj]=0.5F-0.5F*OPUS_COSF((2*OPUS_PI/(_window_sz-1))*xj); |
| } |
| for(xj=0;xj<_window_sz;xj++){ |
| c[xj]=OPUS_COSF((2*OPUS_PI/_window_sz)*xj); |
| } |
| for(xj=0;xj<_window_sz;xj++){ |
| s[xj]=OPUS_SINF((2*OPUS_PI/_window_sz)*xj); |
| } |
| for(xi=0;xi<_nframes;xi++){ |
| int ci; |
| int xk; |
| int bi; |
| for(ci=0;ci<_nchannels;ci++){ |
| for(xk=0;xk<_window_sz;xk++){ |
| x[ci*_window_sz+xk]=window[xk]*_in[(xi*_step+xk)*_nchannels+ci]; |
| } |
| } |
| for(bi=xj=0;bi<_nbands;bi++){ |
| float p[2]={0}; |
| for(;xj<_bands[bi+1];xj++){ |
| for(ci=0;ci<_nchannels;ci++){ |
| float re; |
| float im; |
| int ti; |
| ti=0; |
| re=im=0; |
| for(xk=0;xk<_window_sz;xk++){ |
| re+=c[ti]*x[ci*_window_sz+xk]; |
| im-=s[ti]*x[ci*_window_sz+xk]; |
| ti+=xj; |
| if(ti>=_window_sz)ti-=_window_sz; |
| } |
| re*=_downsample; |
| im*=_downsample; |
| _ps[(xi*ps_sz+xj)*_nchannels+ci]=re*re+im*im+100000; |
| p[ci]+=_ps[(xi*ps_sz+xj)*_nchannels+ci]; |
| } |
| } |
| if(_out){ |
| _out[(xi*_nbands+bi)*_nchannels]=p[0]/(_bands[bi+1]-_bands[bi]); |
| if(_nchannels==2){ |
| _out[(xi*_nbands+bi)*_nchannels+1]=p[1]/(_bands[bi+1]-_bands[bi]); |
| } |
| } |
| } |
| } |
| free(window); |
| } |
| |
| #define NBANDS (21) |
| #define NFREQS (240) |
| |
| /*Bands on which we compute the pseudo-NMR (Bark-derived |
| CELT bands).*/ |
| static const int BANDS[NBANDS+1]={ |
| 0,2,4,6,8,10,12,14,16,20,24,28,32,40,48,56,68,80,96,120,156,200 |
| }; |
| |
| #define TEST_WIN_SIZE (480) |
| #define TEST_WIN_STEP (120) |
| |
| int main(int _argc,const char **_argv){ |
| FILE *fin1; |
| FILE *fin2; |
| float *x; |
| float *y; |
| float *xb; |
| float *X; |
| float *Y; |
| double err; |
| float Q; |
| size_t xlength; |
| size_t ylength; |
| size_t nframes; |
| size_t xi; |
| int ci; |
| int xj; |
| int bi; |
| int nchannels; |
| unsigned rate; |
| int downsample; |
| int ybands; |
| int yfreqs; |
| int max_compare; |
| if(_argc<3||_argc>6){ |
| fprintf(stderr,"Usage: %s [-s] [-r rate2] <file1.sw> <file2.sw>\n", |
| _argv[0]); |
| return EXIT_FAILURE; |
| } |
| nchannels=1; |
| if(strcmp(_argv[1],"-s")==0){ |
| nchannels=2; |
| _argv++; |
| } |
| rate=48000; |
| ybands=NBANDS; |
| yfreqs=NFREQS; |
| downsample=1; |
| if(strcmp(_argv[1],"-r")==0){ |
| rate=atoi(_argv[2]); |
| if(rate!=8000&&rate!=12000&&rate!=16000&&rate!=24000&&rate!=48000){ |
| fprintf(stderr, |
| "Sampling rate must be 8000, 12000, 16000, 24000, or 48000\n"); |
| return EXIT_FAILURE; |
| } |
| downsample=48000/rate; |
| switch(rate){ |
| case 8000:ybands=13;break; |
| case 12000:ybands=15;break; |
| case 16000:ybands=17;break; |
| case 24000:ybands=19;break; |
| } |
| yfreqs=NFREQS/downsample; |
| _argv+=2; |
| } |
| fin1=fopen(_argv[1],"rb"); |
| if(fin1==NULL){ |
| fprintf(stderr,"Error opening '%s'.\n",_argv[1]); |
| return EXIT_FAILURE; |
| } |
| fin2=fopen(_argv[2],"rb"); |
| if(fin2==NULL){ |
| fprintf(stderr,"Error opening '%s'.\n",_argv[2]); |
| fclose(fin1); |
| return EXIT_FAILURE; |
| } |
| /*Read in the data and allocate scratch space.*/ |
| xlength=read_pcm16(&x,fin1,2); |
| if(nchannels==1){ |
| for(xi=0;xi<xlength;xi++)x[xi]=.5*(x[2*xi]+x[2*xi+1]); |
| } |
| fclose(fin1); |
| ylength=read_pcm16(&y,fin2,nchannels); |
| fclose(fin2); |
| if(xlength!=ylength*downsample){ |
| fprintf(stderr,"Sample counts do not match (%lu!=%lu).\n", |
| (unsigned long)xlength,(unsigned long)ylength*downsample); |
| return EXIT_FAILURE; |
| } |
| if(xlength<TEST_WIN_SIZE){ |
| fprintf(stderr,"Insufficient sample data (%lu<%i).\n", |
| (unsigned long)xlength,TEST_WIN_SIZE); |
| return EXIT_FAILURE; |
| } |
| nframes=(xlength-TEST_WIN_SIZE+TEST_WIN_STEP)/TEST_WIN_STEP; |
| xb=(float *)opus_malloc(nframes*NBANDS*nchannels*sizeof(*xb)); |
| X=(float *)opus_malloc(nframes*NFREQS*nchannels*sizeof(*X)); |
| Y=(float *)opus_malloc(nframes*yfreqs*nchannels*sizeof(*Y)); |
| /*Compute the per-band spectral energy of the original signal |
| and the error.*/ |
| band_energy(xb,X,BANDS,NBANDS,x,nchannels,nframes, |
| TEST_WIN_SIZE,TEST_WIN_STEP,1); |
| free(x); |
| band_energy(NULL,Y,BANDS,ybands,y,nchannels,nframes, |
| TEST_WIN_SIZE/downsample,TEST_WIN_STEP/downsample,downsample); |
| free(y); |
| for(xi=0;xi<nframes;xi++){ |
| /*Frequency masking (low to high): 10 dB/Bark slope.*/ |
| for(bi=1;bi<NBANDS;bi++){ |
| for(ci=0;ci<nchannels;ci++){ |
| xb[(xi*NBANDS+bi)*nchannels+ci]+= |
| 0.1F*xb[(xi*NBANDS+bi-1)*nchannels+ci]; |
| } |
| } |
| /*Frequency masking (high to low): 15 dB/Bark slope.*/ |
| for(bi=NBANDS-1;bi-->0;){ |
| for(ci=0;ci<nchannels;ci++){ |
| xb[(xi*NBANDS+bi)*nchannels+ci]+= |
| 0.03F*xb[(xi*NBANDS+bi+1)*nchannels+ci]; |
| } |
| } |
| if(xi>0){ |
| /*Temporal masking: -3 dB/2.5ms slope.*/ |
| for(bi=0;bi<NBANDS;bi++){ |
| for(ci=0;ci<nchannels;ci++){ |
| xb[(xi*NBANDS+bi)*nchannels+ci]+= |
| 0.5F*xb[((xi-1)*NBANDS+bi)*nchannels+ci]; |
| } |
| } |
| } |
| /* Allowing some cross-talk */ |
| if(nchannels==2){ |
| for(bi=0;bi<NBANDS;bi++){ |
| float l,r; |
| l=xb[(xi*NBANDS+bi)*nchannels+0]; |
| r=xb[(xi*NBANDS+bi)*nchannels+1]; |
| xb[(xi*NBANDS+bi)*nchannels+0]+=0.01F*r; |
| xb[(xi*NBANDS+bi)*nchannels+1]+=0.01F*l; |
| } |
| } |
| |
| /* Apply masking */ |
| for(bi=0;bi<ybands;bi++){ |
| for(xj=BANDS[bi];xj<BANDS[bi+1];xj++){ |
| for(ci=0;ci<nchannels;ci++){ |
| X[(xi*NFREQS+xj)*nchannels+ci]+= |
| 0.1F*xb[(xi*NBANDS+bi)*nchannels+ci]; |
| Y[(xi*yfreqs+xj)*nchannels+ci]+= |
| 0.1F*xb[(xi*NBANDS+bi)*nchannels+ci]; |
| } |
| } |
| } |
| } |
| |
| /* Average of consecutive frames to make comparison slightly less sensitive */ |
| for(bi=0;bi<ybands;bi++){ |
| for(xj=BANDS[bi];xj<BANDS[bi+1];xj++){ |
| for(ci=0;ci<nchannels;ci++){ |
| float xtmp; |
| float ytmp; |
| xtmp = X[xj*nchannels+ci]; |
| ytmp = Y[xj*nchannels+ci]; |
| for(xi=1;xi<nframes;xi++){ |
| float xtmp2; |
| float ytmp2; |
| xtmp2 = X[(xi*NFREQS+xj)*nchannels+ci]; |
| ytmp2 = Y[(xi*yfreqs+xj)*nchannels+ci]; |
| X[(xi*NFREQS+xj)*nchannels+ci] += xtmp; |
| Y[(xi*yfreqs+xj)*nchannels+ci] += ytmp; |
| xtmp = xtmp2; |
| ytmp = ytmp2; |
| } |
| } |
| } |
| } |
| |
| /*If working at a lower sampling rate, don't take into account the last |
| 300 Hz to allow for different transition bands. |
| For 12 kHz, we don't skip anything, because the last band already skips |
| 400 Hz.*/ |
| if(rate==48000)max_compare=BANDS[NBANDS]; |
| else if(rate==12000)max_compare=BANDS[ybands]; |
| else max_compare=BANDS[ybands]-3; |
| err=0; |
| for(xi=0;xi<nframes;xi++){ |
| double Ef; |
| Ef=0; |
| for(bi=0;bi<ybands;bi++){ |
| double Eb; |
| Eb=0; |
| for(xj=BANDS[bi];xj<BANDS[bi+1]&&xj<max_compare;xj++){ |
| for(ci=0;ci<nchannels;ci++){ |
| float re; |
| float im; |
| re=Y[(xi*yfreqs+xj)*nchannels+ci]/X[(xi*NFREQS+xj)*nchannels+ci]; |
| im=re-log(re)-1; |
| /*Make comparison less sensitive around the SILK/CELT cross-over to |
| allow for mode freedom in the filters.*/ |
| if(xj>=79&&xj<=81)im*=0.1F; |
| if(xj==80)im*=0.1F; |
| Eb+=im; |
| } |
| } |
| Eb /= (BANDS[bi+1]-BANDS[bi])*nchannels; |
| Ef += Eb*Eb; |
| } |
| /*Using a fixed normalization value means we're willing to accept slightly |
| lower quality for lower sampling rates.*/ |
| Ef/=NBANDS; |
| Ef*=Ef; |
| err+=Ef*Ef; |
| } |
| free(xb); |
| free(X); |
| free(Y); |
| err=pow(err/nframes,1.0/16); |
| Q=100*(1-0.5*log(1+err)/log(1.13)); |
| if(Q<0){ |
| fprintf(stderr,"Test vector FAILS\n"); |
| fprintf(stderr,"Internal weighted error is %f\n",err); |
| return EXIT_FAILURE; |
| } |
| else{ |
| fprintf(stderr,"Test vector PASSES\n"); |
| fprintf(stderr, |
| "Opus quality metric: %.1f %% (internal weighted error is %f)\n",Q,err); |
| return EXIT_SUCCESS; |
| } |
| } |
