Fix pedantic build bustage.
[opus.git] / src / opus_compare.c
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 }