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
38static 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
46static void *opus_malloc(size_t _size){
47  return check_alloc(malloc(_size));
48}
49
50static void *opus_realloc(void *_ptr,size_t _size){
51  return check_alloc(realloc(_ptr,_size));
52}
53
54static 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
88static 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).*/
158static 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
165int 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}
380