Only free state once in kiss-fft failed init path
[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_MIN(_x,_y)      ((_x)<(_y)?(_x):(_y))
9 #define OPUS_MAX(_x,_y)      ((_x)>(_y)?(_x):(_y))
10 #define OPUS_CLAMP(_a,_b,_c) OPUS_MAX(_a,OPUS_MIN(_b,_c))
11 #define OPUS_COSF(_x)        ((float)cos(_x))
12 #define OPUS_SINF(_x)        ((float)sin(_x))
13 #define OPUS_SQRTF(_x)       ((float)sqrt(_x))
14 #define OPUS_LOG10F(_x)      ((float)log10(_x))
15
16 static void *check_alloc(void *_ptr){
17   if(_ptr==NULL){
18     fprintf(stderr,"Out of memory.\n");
19     exit(EXIT_FAILURE);
20   }
21   return _ptr;
22 }
23
24 static void *opus_malloc(size_t _size){
25   return check_alloc(malloc(_size));
26 }
27
28 static void *opus_realloc(void *_ptr,size_t _size){
29   return check_alloc(realloc(_ptr,_size));
30 }
31
32 static size_t read_pcm16(float **_samples,FILE *_fin,
33                          int _nchannels){
34   unsigned char  buf[1024];
35   float         *samples;
36   size_t         nsamples;
37   size_t         csamples;
38   size_t         xi;
39   size_t         nread;
40   samples=NULL;
41   nsamples=csamples=0;
42   for(;;){
43     nread=fread(buf,2*_nchannels,1024/(2*_nchannels),_fin);
44     if(nread<=0)break;
45     if(nsamples+nread>csamples){
46       do csamples=csamples<<1|1;
47       while(nsamples+nread>csamples);
48       samples=(float *)opus_realloc(samples,
49        _nchannels*csamples*sizeof(*samples));
50     }
51     for(xi=0;xi<nread;xi++){
52       int ci;
53       for(ci=0;ci<_nchannels;ci++){
54         int s;
55         s=buf[2*(xi*_nchannels+ci)+1]<<8|buf[2*(xi*_nchannels+ci)];
56         s=((s&0xFFFF)^0x8000)-0x8000;
57         samples[(nsamples+xi)*_nchannels+ci]=s;
58       }
59     }
60     nsamples+=nread;
61   }
62   *_samples=(float *)opus_realloc(samples,
63                      _nchannels*nsamples*sizeof(*samples));
64   return nsamples;
65 }
66
67 static void band_energy(float *_out,const int *_bands,int _nbands,
68  const float *_in,int _nchannels,size_t _nframes,int _window_sz,
69  int _step){
70   float *window;
71   float *x;
72   float *c;
73   float *s;
74   size_t xi;
75   int    xj;
76   window=(float *)opus_malloc((3+_nchannels)*_window_sz
77           *sizeof(*window));
78   c=window+_window_sz;
79   s=c+_window_sz;
80   x=s+_window_sz;
81   for(xj=0;xj<_window_sz;xj++){
82     window[xj]=0.5F-0.5F*OPUS_COSF((2*OPUS_PI/(_window_sz-1))*xj);
83   }
84   for(xj=0;xj<_window_sz;xj++)
85       c[xj]=OPUS_COSF((2*OPUS_PI/_window_sz)*xj);
86   for(xj=0;xj<_window_sz;xj++)
87       s[xj]=OPUS_SINF((2*OPUS_PI/_window_sz)*xj);
88   for(xi=0;xi<_nframes;xi++){
89     int ci;
90     int xk;
91     int bi;
92     for(ci=0;ci<_nchannels;ci++){
93       for(xk=0;xk<_window_sz;xk++){
94         x[ci*_window_sz+xk]=window[xk]
95                            *_in[(xi*_step+xk)*_nchannels+ci];
96       }
97     }
98     for(bi=xj=0;bi<_nbands;bi++){
99       float e2;
100       e2=0;
101       for(;xj<_bands[bi+1];xj++){
102         float p;
103         p=0;
104         for(ci=0;ci<_nchannels;ci++){
105           float re;
106           float im;
107           int   ti;
108           ti=0;
109           re=im=0;
110           for(xk=0;xk<_window_sz;xk++){
111             re+=c[ti]*x[ci*_window_sz+xk];
112             im-=s[ti]*x[ci*_window_sz+xk];
113             ti+=xj;
114             if(ti>=_window_sz)ti-=_window_sz;
115           }
116           p+=OPUS_SQRTF(re*re+im*im);
117         }
118         p*=(1.0F/_nchannels);
119         e2+=p*p;
120       }
121       _out[xi*_nbands+bi]=e2/(_bands[bi+1]-_bands[bi])+1;
122     }
123   }
124   free(window);
125 }
126
127 static int cmp_float(const void *_a,const void *_b){
128   float a;
129   float b;
130   a=*(const float *)_a;
131   b=*(const float *)_b;
132   return (a>b)-(a<b);
133 }
134
135 #define NBANDS (21)
136
137 /*Bands on which we compute the pseudo-NMR (Bark-derived
138   CELT bands).*/
139 static const int BANDS[NBANDS+1]={
140   0,2,4,6,8,10,12,14,16,20,24,28,32,40,48,56,68,80,96,120,156,200
141 };
142
143 /*Per-band NMR threshold.*/
144 static const float NMR_THRESH[NBANDS]={
145 85113.804F,72443.596F,61659.5F,  52480.746F,44668.359F,38018.940F,
146 32359.366F,27542.287F,23442.288F,19952.623F,16982.437F,14454.398F,
147 12302.688F,10471.285F, 8912.5094F,7585.7758F,6456.5423F,5495.4087F,
148  4677.3514F,3981.0717F,3388.4416F
149 };
150
151 /*Noise floor.*/
152 static const float NOISE_FLOOR[NBANDS]={
153 8.7096359F,7.5857758F,6.6069345F,5.7543994F,5.0118723F,4.3651583F,
154 3.8018940F,3.3113112F,2.8840315F,2.5118864F,2.1877616F,1.9054607F,
155 1.6595869F,1.4454398F,1.2589254F,1.0964782F,0.95499259F,0.83176377F,
156 0.72443596F,0.63095734F,0.54954087F
157 };
158
159 #define TEST_WIN_SIZE (480)
160 #define TEST_WIN_STEP (TEST_WIN_SIZE>>1)
161
162 int main(int _argc,const char **_argv){
163   FILE   *fin1;
164   FILE   *fin2;
165   float  *x;
166   float  *y;
167   float  *xb;
168   float  *eb;
169   float  *nmr;
170   float   thresh;
171   float   mismatch;
172   float   err;
173   float   nmr_sum;
174   size_t  weight;
175   size_t  xlength;
176   size_t  ylength;
177   size_t  nframes;
178   size_t  xi;
179   int     bi;
180   int     nchannels;
181   if(_argc<3||_argc>4){
182     fprintf(stderr,"Usage: %s [-s] <file1.sw> <file2.sw>\n",
183             _argv[0]);
184     return EXIT_FAILURE;
185   }
186   nchannels=1;
187   if(strcmp(_argv[1],"-s")==0)nchannels=2;
188   fin1=fopen(_argv[nchannels],"rb");
189   if(fin1==NULL){
190     fprintf(stderr,"Error opening '%s'.\n",_argv[nchannels]);
191     return EXIT_FAILURE;
192   }
193   fin2=fopen(_argv[nchannels+1],"rb");
194   if(fin2==NULL){
195     fprintf(stderr,"Error opening '%s'.\n",_argv[nchannels+1]);
196     fclose(fin1);
197     return EXIT_FAILURE;
198   }
199   /*Read in the data and allocate scratch space.*/
200   xlength=read_pcm16(&x,fin1,nchannels);
201   fclose(fin1);
202   ylength=read_pcm16(&y,fin2,nchannels);
203   fclose(fin2);
204   if(xlength!=ylength){
205     fprintf(stderr,"Sample counts do not match (%lu!=%lu).\n",
206      (unsigned long)xlength,(unsigned long)ylength);
207     return EXIT_FAILURE;
208   }
209   if(xlength<TEST_WIN_SIZE){
210     fprintf(stderr,"Insufficient sample data (%lu<%i).\n",
211      (unsigned long)xlength,TEST_WIN_SIZE);
212     return EXIT_FAILURE;
213   }
214   nframes=(xlength-TEST_WIN_SIZE+TEST_WIN_STEP)/TEST_WIN_STEP;
215   xb=(float *)opus_malloc(nframes*NBANDS*sizeof(*xb));
216   eb=(float *)opus_malloc(nframes*NBANDS*sizeof(*eb));
217   nmr=(float *)opus_malloc(nframes*NBANDS*sizeof(*nmr));
218   /*Compute the error signal.*/
219   for(xi=0;xi<xlength*nchannels;xi++){
220     err=x[xi]-y[xi];
221     y[xi]=err-OPUS_CLAMP(-1,err,1);
222   }
223   /*Compute the per-band spectral energy of the original signal
224     and the error.*/
225   band_energy(xb,BANDS,NBANDS,x,nchannels,nframes,
226           TEST_WIN_SIZE,TEST_WIN_STEP);
227   free(x);
228   band_energy(eb,BANDS,NBANDS,y,nchannels,nframes,
229           TEST_WIN_SIZE,TEST_WIN_STEP);
230   free(y);
231   nmr_sum=0;
232   for(xi=0;xi<nframes;xi++){
233     /*Frequency masking (low to high): 10 dB/Bark slope.*/
234     for(bi=1;bi<NBANDS;bi++)
235         xb[xi*NBANDS+bi]+=0.1F*xb[xi*NBANDS+bi-1];
236     /*Frequency masking (high to low): 15 dB/Bark slope.*/
237     for(bi=NBANDS-1;bi-->0;)
238         xb[xi*NBANDS+bi]+=0.03F*xb[xi*NBANDS+bi+1];
239     if(xi>0){
240       /*Temporal masking: 5 dB/5ms slope.*/
241       for(bi=0;bi<NBANDS;bi++)
242           xb[xi*NBANDS+bi]+=0.3F*xb[(xi-1)*NBANDS+bi];
243     }
244     /*Compute NMR.*/
245     for(bi=0;bi<NBANDS;bi++){
246       nmr[xi*NBANDS+bi]=xb[xi*NBANDS+bi]/eb[xi*NBANDS+bi];
247       nmr_sum+=10*OPUS_LOG10F(nmr[xi*NBANDS+bi]);
248     }
249   }
250   /*Find the 90th percentile of the errors.*/
251   memcpy(xb,eb,nframes*NBANDS*sizeof(*xb));
252   qsort(xb,nframes*NBANDS,sizeof(*xb),cmp_float);
253   thresh=xb[(9*nframes*NBANDS+5)/10];
254   free(xb);
255   /*Compute the mismatch.*/
256   mismatch=0;
257   weight=0;
258   for(xi=0;xi<nframes;xi++){
259     for(bi=0;bi<NBANDS;bi++){
260       if(eb[xi*NBANDS+bi]>thresh){
261         mismatch+=NMR_THRESH[bi]/nmr[xi*NBANDS+bi];
262         weight++;
263       }
264     }
265   }
266   free(nmr);
267   free(eb);
268   printf("Average pseudo-NMR: %3.2f dB\n",nmr_sum/(nframes*NBANDS));
269   if(weight<=0){
270     err=-100;
271     printf("Mismatch level: below noise floor\n");
272   }
273   else{
274     err=10*OPUS_LOG10F(mismatch/weight);
275     printf("Weighted mismatch: %3.2f dB\n",err);
276   }
277   printf("\n");
278   if(err<0){
279     printf("**Decoder PASSES test (mismatch < 0 dB)\n");
280     return EXIT_SUCCESS;
281   }
282   printf("**Decoder FAILS test (mismatch >= 0 dB)\n");
283   return EXIT_FAILURE;
284 }