Fix MSVC linker warnings
[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 (120)
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   double    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: -3 dB/2.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.5F*xb[((xi-1)*NBANDS+bi)*nchannels+ci];
254         }
255       }
256     }
257     /* Allowing some cross-talk */
258     if(nchannels==2){
259       for(bi=0;bi<NBANDS;bi++){
260         float l,r;
261         l=xb[(xi*NBANDS+bi)*nchannels+0];
262         r=xb[(xi*NBANDS+bi)*nchannels+1];
263         xb[(xi*NBANDS+bi)*nchannels+0]+=0.01F*r;
264         xb[(xi*NBANDS+bi)*nchannels+1]+=0.01F*l;
265       }
266     }
267
268     /* Apply masking */
269     for(bi=0;bi<ybands;bi++){
270       for(xj=BANDS[bi];xj<BANDS[bi+1];xj++){
271         for(ci=0;ci<nchannels;ci++){
272           X[(xi*NFREQS+xj)*nchannels+ci]+=
273            0.1F*xb[(xi*NBANDS+bi)*nchannels+ci];
274           Y[(xi*yfreqs+xj)*nchannels+ci]+=
275            0.1F*xb[(xi*NBANDS+bi)*nchannels+ci];
276         }
277       }
278     }
279   }
280
281   /* Average of consecutive frames to make comparison slightly less sensitive */
282   for(bi=0;bi<ybands;bi++){
283     for(xj=BANDS[bi];xj<BANDS[bi+1];xj++){
284       for(ci=0;ci<nchannels;ci++){
285          float xtmp;
286          float ytmp;
287          xtmp = X[xj*nchannels+ci];
288          ytmp = Y[xj*nchannels+ci];
289          for(xi=1;xi<nframes;xi++){
290            float xtmp2;
291            float ytmp2;
292            xtmp2 = X[(xi*NFREQS+xj)*nchannels+ci];
293            ytmp2 = Y[(xi*yfreqs+xj)*nchannels+ci];
294            X[(xi*NFREQS+xj)*nchannels+ci] += xtmp;
295            Y[(xi*yfreqs+xj)*nchannels+ci] += ytmp;
296            xtmp = xtmp2;
297            ytmp = ytmp2;
298          }
299       }
300     }
301   }
302
303   /*If working at a lower sampling rate, don't take into account the last
304      300 Hz to allow for different transition bands.
305     For 12 kHz, we don't skip anything, because the last band already skips
306      400 Hz.*/
307   if(rate==48000)max_compare=BANDS[NBANDS];
308   else if(rate==12000)max_compare=BANDS[ybands];
309   else max_compare=BANDS[ybands]-3;
310   err=0;
311   for(xi=0;xi<nframes;xi++){
312     double Ef;
313     Ef=0;
314     for(bi=0;bi<ybands;bi++){
315       double Eb;
316       Eb=0;
317       for(xj=BANDS[bi];xj<BANDS[bi+1]&&xj<max_compare;xj++){
318         for(ci=0;ci<nchannels;ci++){
319           float re;
320           float im;
321           re=Y[(xi*yfreqs+xj)*nchannels+ci]/X[(xi*NFREQS+xj)*nchannels+ci];
322           im=re-log(re)-1;
323           /*Make comparison less sensitive around the SILK/CELT cross-over to
324             allow for mode freedom in the filters.*/
325           if(xj>=79&&xj<=81)im*=0.1F;
326           if(xj==80)im*=0.1F;
327           Eb+=im;
328         }
329       }
330       Eb /= (BANDS[bi+1]-BANDS[bi])*nchannels;
331       Ef += Eb*Eb;
332     }
333     /*Using a fixed normalization value means we're willing to accept slightly
334        lower quality for lower sampling rates.*/
335     Ef/=NBANDS;
336     Ef*=Ef;
337     err+=Ef*Ef;
338   }
339   err=pow(err/nframes,1.0/16);
340   Q=100*(1-0.5*log(1+err)/log(1.13));
341   if(Q<0){
342     fprintf(stderr,"Test vector FAILS\n");
343     fprintf(stderr,"Internal weighted error is %f\n",err);
344     return EXIT_FAILURE;
345   }
346   else{
347     fprintf(stderr,"Test vector PASSES\n");
348     fprintf(stderr,
349      "Opus quality metric: %.1f %% (internal weighted error is %f)\n",Q,err);
350     return EXIT_SUCCESS;
351   }
352 }