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