382 lines
11 KiB
C
382 lines
11 KiB
C
/* Copyright (c) 2011-2012 Xiph.Org Foundation, Mozilla Corporation
|
|
Written by Jean-Marc Valin and Timothy B. Terriberry */
|
|
/*
|
|
Redistribution and use in source and binary forms, with or without
|
|
modification, are permitted provided that the following conditions
|
|
are met:
|
|
|
|
- Redistributions of source code must retain the above copyright
|
|
notice, this list of conditions and the following disclaimer.
|
|
|
|
- Redistributions in binary form must reproduce the above copyright
|
|
notice, this list of conditions and the following disclaimer in the
|
|
documentation and/or other materials provided with the distribution.
|
|
|
|
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
|
|
``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
|
|
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
|
|
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
|
|
OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
|
|
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
|
|
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
|
|
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
|
|
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
|
|
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
|
|
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
|
*/
|
|
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <math.h>
|
|
#include <string.h>
|
|
|
|
#define OPUS_PI (3.14159265F)
|
|
|
|
#define OPUS_COSF(_x) ((float)cos(_x))
|
|
#define OPUS_SINF(_x) ((float)sin(_x))
|
|
|
|
static void *check_alloc(void *_ptr){
|
|
if(_ptr==NULL){
|
|
fprintf(stderr,"Out of memory.\n");
|
|
exit(EXIT_FAILURE);
|
|
}
|
|
return _ptr;
|
|
}
|
|
|
|
static void *opus_malloc(size_t _size){
|
|
return check_alloc(malloc(_size));
|
|
}
|
|
|
|
static void *opus_realloc(void *_ptr,size_t _size){
|
|
return check_alloc(realloc(_ptr,_size));
|
|
}
|
|
|
|
static size_t read_pcm16(float **_samples,FILE *_fin,int _nchannels){
|
|
unsigned char buf[1024];
|
|
float *samples;
|
|
size_t nsamples;
|
|
size_t csamples;
|
|
size_t xi;
|
|
size_t nread;
|
|
samples=NULL;
|
|
nsamples=csamples=0;
|
|
for(;;){
|
|
nread=fread(buf,2*_nchannels,1024/(2*_nchannels),_fin);
|
|
if(nread<=0)break;
|
|
if(nsamples+nread>csamples){
|
|
do csamples=csamples<<1|1;
|
|
while(nsamples+nread>csamples);
|
|
samples=(float *)opus_realloc(samples,
|
|
_nchannels*csamples*sizeof(*samples));
|
|
}
|
|
for(xi=0;xi<nread;xi++){
|
|
int ci;
|
|
for(ci=0;ci<_nchannels;ci++){
|
|
int s;
|
|
s=buf[2*(xi*_nchannels+ci)+1]<<8|buf[2*(xi*_nchannels+ci)];
|
|
s=((s&0xFFFF)^0x8000)-0x8000;
|
|
samples[(nsamples+xi)*_nchannels+ci]=s;
|
|
}
|
|
}
|
|
nsamples+=nread;
|
|
}
|
|
*_samples=(float *)opus_realloc(samples,
|
|
_nchannels*nsamples*sizeof(*samples));
|
|
return nsamples;
|
|
}
|
|
|
|
static void band_energy(float *_out,float *_ps,const int *_bands,int _nbands,
|
|
const float *_in,int _nchannels,size_t _nframes,int _window_sz,
|
|
int _step,int _downsample){
|
|
float *window;
|
|
float *x;
|
|
float *c;
|
|
float *s;
|
|
size_t xi;
|
|
int xj;
|
|
int ps_sz;
|
|
window=(float *)opus_malloc((3+_nchannels)*_window_sz*sizeof(*window));
|
|
c=window+_window_sz;
|
|
s=c+_window_sz;
|
|
x=s+_window_sz;
|
|
ps_sz=_window_sz/2;
|
|
for(xj=0;xj<_window_sz;xj++){
|
|
window[xj]=0.5F-0.5F*OPUS_COSF((2*OPUS_PI/(_window_sz-1))*xj);
|
|
}
|
|
for(xj=0;xj<_window_sz;xj++){
|
|
c[xj]=OPUS_COSF((2*OPUS_PI/_window_sz)*xj);
|
|
}
|
|
for(xj=0;xj<_window_sz;xj++){
|
|
s[xj]=OPUS_SINF((2*OPUS_PI/_window_sz)*xj);
|
|
}
|
|
for(xi=0;xi<_nframes;xi++){
|
|
int ci;
|
|
int xk;
|
|
int bi;
|
|
for(ci=0;ci<_nchannels;ci++){
|
|
for(xk=0;xk<_window_sz;xk++){
|
|
x[ci*_window_sz+xk]=window[xk]*_in[(xi*_step+xk)*_nchannels+ci];
|
|
}
|
|
}
|
|
for(bi=xj=0;bi<_nbands;bi++){
|
|
float p[2]={0};
|
|
for(;xj<_bands[bi+1];xj++){
|
|
for(ci=0;ci<_nchannels;ci++){
|
|
float re;
|
|
float im;
|
|
int ti;
|
|
ti=0;
|
|
re=im=0;
|
|
for(xk=0;xk<_window_sz;xk++){
|
|
re+=c[ti]*x[ci*_window_sz+xk];
|
|
im-=s[ti]*x[ci*_window_sz+xk];
|
|
ti+=xj;
|
|
if(ti>=_window_sz)ti-=_window_sz;
|
|
}
|
|
re*=_downsample;
|
|
im*=_downsample;
|
|
_ps[(xi*ps_sz+xj)*_nchannels+ci]=re*re+im*im+100000;
|
|
p[ci]+=_ps[(xi*ps_sz+xj)*_nchannels+ci];
|
|
}
|
|
}
|
|
if(_out){
|
|
_out[(xi*_nbands+bi)*_nchannels]=p[0]/(_bands[bi+1]-_bands[bi]);
|
|
if(_nchannels==2){
|
|
_out[(xi*_nbands+bi)*_nchannels+1]=p[1]/(_bands[bi+1]-_bands[bi]);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
free(window);
|
|
}
|
|
|
|
#define NBANDS (21)
|
|
#define NFREQS (240)
|
|
|
|
/*Bands on which we compute the pseudo-NMR (Bark-derived
|
|
CELT bands).*/
|
|
static const int BANDS[NBANDS+1]={
|
|
0,2,4,6,8,10,12,14,16,20,24,28,32,40,48,56,68,80,96,120,156,200
|
|
};
|
|
|
|
#define TEST_WIN_SIZE (480)
|
|
#define TEST_WIN_STEP (120)
|
|
|
|
int main(int _argc,const char **_argv){
|
|
FILE *fin1;
|
|
FILE *fin2;
|
|
float *x;
|
|
float *y;
|
|
float *xb;
|
|
float *X;
|
|
float *Y;
|
|
double err;
|
|
float Q;
|
|
size_t xlength;
|
|
size_t ylength;
|
|
size_t nframes;
|
|
size_t xi;
|
|
int ci;
|
|
int xj;
|
|
int bi;
|
|
int nchannels;
|
|
unsigned rate;
|
|
int downsample;
|
|
int ybands;
|
|
int yfreqs;
|
|
int max_compare;
|
|
if(_argc<3||_argc>6){
|
|
fprintf(stderr,"Usage: %s [-s] [-r rate2] <file1.sw> <file2.sw>\n",
|
|
_argv[0]);
|
|
return EXIT_FAILURE;
|
|
}
|
|
nchannels=1;
|
|
if(strcmp(_argv[1],"-s")==0){
|
|
nchannels=2;
|
|
_argv++;
|
|
}
|
|
rate=48000;
|
|
ybands=NBANDS;
|
|
yfreqs=NFREQS;
|
|
downsample=1;
|
|
if(strcmp(_argv[1],"-r")==0){
|
|
rate=atoi(_argv[2]);
|
|
if(rate!=8000&&rate!=12000&&rate!=16000&&rate!=24000&&rate!=48000){
|
|
fprintf(stderr,
|
|
"Sampling rate must be 8000, 12000, 16000, 24000, or 48000\n");
|
|
return EXIT_FAILURE;
|
|
}
|
|
downsample=48000/rate;
|
|
switch(rate){
|
|
case 8000:ybands=13;break;
|
|
case 12000:ybands=15;break;
|
|
case 16000:ybands=17;break;
|
|
case 24000:ybands=19;break;
|
|
}
|
|
yfreqs=NFREQS/downsample;
|
|
_argv+=2;
|
|
}
|
|
fin1=fopen(_argv[1],"rb");
|
|
if(fin1==NULL){
|
|
fprintf(stderr,"Error opening '%s'.\n",_argv[1]);
|
|
return EXIT_FAILURE;
|
|
}
|
|
fin2=fopen(_argv[2],"rb");
|
|
if(fin2==NULL){
|
|
fprintf(stderr,"Error opening '%s'.\n",_argv[2]);
|
|
fclose(fin1);
|
|
return EXIT_FAILURE;
|
|
}
|
|
/*Read in the data and allocate scratch space.*/
|
|
xlength=read_pcm16(&x,fin1,2);
|
|
if(nchannels==1){
|
|
for(xi=0;xi<xlength;xi++)x[xi]=.5*(x[2*xi]+x[2*xi+1]);
|
|
}
|
|
fclose(fin1);
|
|
ylength=read_pcm16(&y,fin2,nchannels);
|
|
fclose(fin2);
|
|
if(xlength!=ylength*downsample){
|
|
fprintf(stderr,"Sample counts do not match (%lu!=%lu).\n",
|
|
(unsigned long)xlength,(unsigned long)ylength*downsample);
|
|
return EXIT_FAILURE;
|
|
}
|
|
if(xlength<TEST_WIN_SIZE){
|
|
fprintf(stderr,"Insufficient sample data (%lu<%i).\n",
|
|
(unsigned long)xlength,TEST_WIN_SIZE);
|
|
return EXIT_FAILURE;
|
|
}
|
|
nframes=(xlength-TEST_WIN_SIZE+TEST_WIN_STEP)/TEST_WIN_STEP;
|
|
xb=(float *)opus_malloc(nframes*NBANDS*nchannels*sizeof(*xb));
|
|
X=(float *)opus_malloc(nframes*NFREQS*nchannels*sizeof(*X));
|
|
Y=(float *)opus_malloc(nframes*yfreqs*nchannels*sizeof(*Y));
|
|
/*Compute the per-band spectral energy of the original signal
|
|
and the error.*/
|
|
band_energy(xb,X,BANDS,NBANDS,x,nchannels,nframes,
|
|
TEST_WIN_SIZE,TEST_WIN_STEP,1);
|
|
free(x);
|
|
band_energy(NULL,Y,BANDS,ybands,y,nchannels,nframes,
|
|
TEST_WIN_SIZE/downsample,TEST_WIN_STEP/downsample,downsample);
|
|
free(y);
|
|
for(xi=0;xi<nframes;xi++){
|
|
/*Frequency masking (low to high): 10 dB/Bark slope.*/
|
|
for(bi=1;bi<NBANDS;bi++){
|
|
for(ci=0;ci<nchannels;ci++){
|
|
xb[(xi*NBANDS+bi)*nchannels+ci]+=
|
|
0.1F*xb[(xi*NBANDS+bi-1)*nchannels+ci];
|
|
}
|
|
}
|
|
/*Frequency masking (high to low): 15 dB/Bark slope.*/
|
|
for(bi=NBANDS-1;bi-->0;){
|
|
for(ci=0;ci<nchannels;ci++){
|
|
xb[(xi*NBANDS+bi)*nchannels+ci]+=
|
|
0.03F*xb[(xi*NBANDS+bi+1)*nchannels+ci];
|
|
}
|
|
}
|
|
if(xi>0){
|
|
/*Temporal masking: -3 dB/2.5ms slope.*/
|
|
for(bi=0;bi<NBANDS;bi++){
|
|
for(ci=0;ci<nchannels;ci++){
|
|
xb[(xi*NBANDS+bi)*nchannels+ci]+=
|
|
0.5F*xb[((xi-1)*NBANDS+bi)*nchannels+ci];
|
|
}
|
|
}
|
|
}
|
|
/* Allowing some cross-talk */
|
|
if(nchannels==2){
|
|
for(bi=0;bi<NBANDS;bi++){
|
|
float l,r;
|
|
l=xb[(xi*NBANDS+bi)*nchannels+0];
|
|
r=xb[(xi*NBANDS+bi)*nchannels+1];
|
|
xb[(xi*NBANDS+bi)*nchannels+0]+=0.01F*r;
|
|
xb[(xi*NBANDS+bi)*nchannels+1]+=0.01F*l;
|
|
}
|
|
}
|
|
|
|
/* Apply masking */
|
|
for(bi=0;bi<ybands;bi++){
|
|
for(xj=BANDS[bi];xj<BANDS[bi+1];xj++){
|
|
for(ci=0;ci<nchannels;ci++){
|
|
X[(xi*NFREQS+xj)*nchannels+ci]+=
|
|
0.1F*xb[(xi*NBANDS+bi)*nchannels+ci];
|
|
Y[(xi*yfreqs+xj)*nchannels+ci]+=
|
|
0.1F*xb[(xi*NBANDS+bi)*nchannels+ci];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/* Average of consecutive frames to make comparison slightly less sensitive */
|
|
for(bi=0;bi<ybands;bi++){
|
|
for(xj=BANDS[bi];xj<BANDS[bi+1];xj++){
|
|
for(ci=0;ci<nchannels;ci++){
|
|
float xtmp;
|
|
float ytmp;
|
|
xtmp = X[xj*nchannels+ci];
|
|
ytmp = Y[xj*nchannels+ci];
|
|
for(xi=1;xi<nframes;xi++){
|
|
float xtmp2;
|
|
float ytmp2;
|
|
xtmp2 = X[(xi*NFREQS+xj)*nchannels+ci];
|
|
ytmp2 = Y[(xi*yfreqs+xj)*nchannels+ci];
|
|
X[(xi*NFREQS+xj)*nchannels+ci] += xtmp;
|
|
Y[(xi*yfreqs+xj)*nchannels+ci] += ytmp;
|
|
xtmp = xtmp2;
|
|
ytmp = ytmp2;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/*If working at a lower sampling rate, don't take into account the last
|
|
300 Hz to allow for different transition bands.
|
|
For 12 kHz, we don't skip anything, because the last band already skips
|
|
400 Hz.*/
|
|
if(rate==48000)max_compare=BANDS[NBANDS];
|
|
else if(rate==12000)max_compare=BANDS[ybands];
|
|
else max_compare=BANDS[ybands]-3;
|
|
err=0;
|
|
for(xi=0;xi<nframes;xi++){
|
|
double Ef;
|
|
Ef=0;
|
|
for(bi=0;bi<ybands;bi++){
|
|
double Eb;
|
|
Eb=0;
|
|
for(xj=BANDS[bi];xj<BANDS[bi+1]&&xj<max_compare;xj++){
|
|
for(ci=0;ci<nchannels;ci++){
|
|
float re;
|
|
float im;
|
|
re=Y[(xi*yfreqs+xj)*nchannels+ci]/X[(xi*NFREQS+xj)*nchannels+ci];
|
|
im=re-log(re)-1;
|
|
/*Make comparison less sensitive around the SILK/CELT cross-over to
|
|
allow for mode freedom in the filters.*/
|
|
if(xj>=79&&xj<=81)im*=0.1F;
|
|
if(xj==80)im*=0.1F;
|
|
Eb+=im;
|
|
}
|
|
}
|
|
Eb /= (BANDS[bi+1]-BANDS[bi])*nchannels;
|
|
Ef += Eb*Eb;
|
|
}
|
|
/*Using a fixed normalization value means we're willing to accept slightly
|
|
lower quality for lower sampling rates.*/
|
|
Ef/=NBANDS;
|
|
Ef*=Ef;
|
|
err+=Ef*Ef;
|
|
}
|
|
free(xb);
|
|
free(X);
|
|
free(Y);
|
|
err=pow(err/nframes,1.0/16);
|
|
Q=100*(1-0.5*log(1+err)/log(1.13));
|
|
if(Q<0){
|
|
fprintf(stderr,"Test vector FAILS\n");
|
|
fprintf(stderr,"Internal weighted error is %f\n",err);
|
|
return EXIT_FAILURE;
|
|
}
|
|
else{
|
|
fprintf(stderr,"Test vector PASSES\n");
|
|
fprintf(stderr,
|
|
"Opus quality metric: %.1f %% (internal weighted error is %f)\n",Q,err);
|
|
return EXIT_SUCCESS;
|
|
}
|
|
}
|