lib/psytune.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
   5  * THE GNU LESSER/LIBRARY PUBLIC LICENSE, WHICH IS INCLUDED WITH    *
   6  * THIS SOURCE. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.        *
   7  *                                                                  *
   8  * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2000             *
   9  * by Monty <monty@xiph.org> and the XIPHOPHORUS Company            *
  10  * http://www.xiph.org/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14  function: simple utility that runs audio through the psychoacoustics
  15            without encoding
  16  last mod: $Id: psytune.c,v 1.12 2001/01/22 01:38:25 xiphmont Exp $
  17
  18  ********************************************************************/
  19
  20 #include <stdio.h>
  21 #include <stdlib.h>
  22 #include <string.h>
  23 #include <math.h>
  24
  25 #include "vorbis/codec.h"
  26 #include "os.h"
  27 #include "psy.h"
  28 #include "mdct.h"
  29 #include "smallft.h"
  30 #include "window.h"
  31 #include "scales.h"
  32 #include "lpc.h"
  33 #include "lsp.h"
  34
  35 static vorbis_info_psy _psy_set0={
  36   1,/*athp*/
  37   0,/*decayp*/
  38
  39   -100.f,
  40   -140.f,
  41
  42   8,
  43
  44   /*     0  1  2   3   4   5   6   7   8   9  10  11  12  13  14  15   16   */
  45   /* x: 63 88 125 175 250 350 500 700 1k 1.4k 2k 2.8k 4k 5.6k 8k 11.5k 16k Hz */
  46   /* y: 0 10 20 30 40 50 60 70 80 90 100 dB */
  47
  48    1,/* tonemaskp */
  49   /*  0   10   20   30   40   50   60   70   80   90   100 */
  50   {
  51    {-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f}, /*63*/
  52    {-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f}, /*88*/
  53    {-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f}, /*125*/
  54    //   {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-99.,-99.,-99.,-100.}, /*175*/
  55    //   {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-99.,-99.,-99.,-100.}, /*250*/
  56    //   {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-99.,-99.,-99.,-100.}, /*350*/
  57    //   {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-99.,-99.,-99.,-100.}, /*500*/
  58    //   {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-99.,-99.,-99.,-100.}, /*700*/
  59
  60    //   {-30.,-35.,-35.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*63*/
  61    //   {-30.,-35.,-35.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*88*/
  62    //   {-30.,-35.,-35.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*125*/
  63    {-30.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*175*/
  64    {-30.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*250*/
  65    {-30.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*350*/
  66    {-30.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*500*/
  67    {-30.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*700*/
  68
  69    {-30.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*1000*/
  70    {-30.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*1400*/
  71    {-40.f,-40.f,-40.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*2000*/
  72    {-40.f,-40.f,-40.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*2800*/
  73    {-30.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*4000*/
  74    {-30.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*5600*/
  75
  76    {-30.f,-30.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*8000*/
  77    {-30.f,-30.f,-30.f,-40.f,-40.f,-45.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*11500*/
  78    {-30.f,-30.f,-30.f,-35.f,-35.f,-45.f,-50.f,-60.f,-70.f,-80.f,-90.f}, /*16000*/
  79   },
  80
  81   1,/* peakattp */
  82   {{-14.f,-16.f,-18.f,-19.f,-20.f,-21.f,-22.f,-22.f,-24.f,-24.f,-24.f},/*63*/
  83    {-14.f,-16.f,-18.f,-19.f,-20.f,-21.f,-22.f,-22.f,-24.f,-24.f,-24.f},/*88*/
  84    {-14.f,-16.f,-18.f,-19.f,-20.f,-21.f,-22.f,-22.f,-24.f,-24.f,-24.f},/*125*/
  85    {-10.f,-10.f,-10.f,-10.f,-16.f,-16.f,-18.f,-20.f,-24.f,-24.f,-24.f},/*175*/
  86    {-10.f,-10.f,-10.f,-10.f,-16.f,-16.f,-18.f,-20.f,-24.f,-24.f,-24.f},/*250*/
  87    {-10.f,-10.f,-10.f,-10.f,-16.f,-16.f,-18.f,-20.f,-22.f,-24.f,-24.f},/*350*/
  88    {-10.f,-10.f,-10.f,-10.f,-16.f,-16.f,-18.f,-20.f,-22.f,-24.f,-24.f},/*500*/
  89    {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*700*/
  90    {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*1000*/
  91    {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*1400*/
  92    {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*2000*/
  93    {-10.f,-10.f,-10.f,-12.f,-16.f,-16.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*2400*/
  94    {-10.f,-10.f,-10.f,-12.f,-16.f,-16.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*4000*/
  95    {-10.f,-10.f,-10.f,-12.f,-12.f,-14.f,-16.f,-18.f,-22.f,-24.f,-24.f},/*5600*/
  96    {-10.f,-10.f,-10.f,-10.f,-10.f,-14.f,-16.f,-18.f,-22.f,-24.f,-24.f},/*8000*/
  97    {-10.f,-10.f,-10.f,-10.f,-10.f,-14.f,-16.f,-18.f,-22.f,-24.f,-24.f},/*11500*/
  98    {-10.f,-10.f,-10.f,-10.f,-10.f,-12.f,-16.f,-18.f,-22.f,-24.f,-24.f},/*16000*/
  99   },
 100
 101   1,/*noisemaskp */
 102   -0.f,  /* suppress any noise curve over maxspec+n */
 103   .5f,   /* low window */
 104   .5f,   /* high window */
 105   25,
 106   25,
 107   {.000f, /*63*/
 108    .000f, /*88*/
 109    .000f, /*125*/
 110    .000f, /*175*/
 111    .000f, /*250*/
 112    .000f, /*350*/
 113    .000f, /*500*/
 114    .200f, /*700*/
 115    .300f, /*1000*/
 116    .400f, /*1400*/
 117    .400f, /*2000*/
 118    .400f, /*2800*/
 119    .700f, /*4000*/
 120    .900f, /*5600*/
 121    .900f, /*8000*/
 122    .900f, /*11500*/
 123    .900f, /*16000*/
 124   },
 125
 126   105.f,  /* even decade + 5 is important; saves an rint() later in a
 127             tight loop) */
 128   -28.,
 129
 130   -0.f, -.004f   /* attack/decay control */
 131 };
 132
 133 static int noisy=0;
 134 void analysis(char *base,int i,float *v,int n,int bark,int dB){
 135   if(noisy){
 136     int j;
 137     FILE *of;
 138     char buffer[80];
 139     sprintf(buffer,"%s_%d.m",base,i);
 140     of=fopen(buffer,"w");
 141
 142     for(j=0;j<n;j++){
 143       if(dB && v[j]==0)
 144           fprintf(of,"\n\n");
 145       else{
 146         if(bark)
 147           fprintf(of,"%g ",toBARK(22050.f*j/n));
 148         else
 149           fprintf(of,"%g ",(float)j);
 150
 151         if(dB){
 152           fprintf(of,"%g\n",todB(fabs(v[j])));
 153         }else{
 154           fprintf(of,"%g\n",v[j]);
 155         }
 156       }
 157     }
 158     fclose(of);
 159   }
 160 }
 161
 162 typedef struct {
 163   long n;
 164   int ln;
 165   int  m;
 166   int *linearmap;
 167
 168   void *vi;
 169   lpc_lookup lpclook;
 170   float *lsp_look;
 171
 172 } vorbis_look_floor0;
 173
 174 extern float _curve_to_lpc(float *curve,float *lpc,
 175                            vorbis_look_floor0 *l);
 176
 177 long frameno=0;
 178
 179 /* hacked from floor0.c */
 180 static void floorinit(vorbis_look_floor0 *look,int n,int m,int ln){
 181   int j;
 182   float scale;
 183   look->m=m;
 184   look->n=n;
 185   look->ln=ln;
 186   lpc_init(&look->lpclook,look->ln,look->m);
 187
 188   scale=look->ln/toBARK(22050.f);
 189
 190   look->linearmap=_ogg_malloc(look->n*sizeof(int));
 191   for(j=0;j<look->n;j++){
 192     int val=floor( toBARK(22050.f/n*j) *scale);
 193     if(val>look->ln)val=look->ln;
 194     look->linearmap[j]=val;
 195   }
 196 }
 197
 198 int main(int argc,char *argv[]){
 199   int eos=0;
 200   float nonz=0.f;
 201   float acc=0.f;
 202   float tot=0.f;
 203   float ampmax=-9999,newmax;
 204
 205   int framesize=2048;
 206   int order=30;
 207   int map=256;
 208   float ampmax_att_per_sec=-10.;
 209
 210   float *pcm[2],*out[2],*window,*lpc,*flr,*mask;
 211   signed char *buffer,*buffer2;
 212   mdct_lookup m_look;
 213   drft_lookup f_look;
 214   vorbis_look_psy p_look;
 215   long i,j,k;
 216
 217   vorbis_look_floor0 floorlook;
 218
 219   int ath=0;
 220   int decayp=0;
 221
 222   argv++;
 223   while(*argv){
 224     if(*argv[0]=='-'){
 225       /* option */
 226       if(argv[0][1]=='v'){
 227         noisy=0;
 228       }
 229       if(argv[0][1]=='o'){
 230         order=atoi(argv[0]+2);
 231       }
 232       if(argv[0][1]=='m'){
 233         map=atoi(argv[0]+2);
 234       }
 235     }else
 236       if(*argv[0]=='+'){
 237         /* option */
 238         if(argv[0][1]=='v'){
 239           noisy=1;
 240         }
 241         if(argv[0][1]=='o'){
 242           order=atoi(argv[0]+2);
 243         }
 244         if(argv[0][1]=='m'){
 245           map=atoi(argv[0]+2);
 246         }
 247       }else
 248         framesize=atoi(argv[0]);
 249     argv++;
 250   }
 251
 252   mask=_ogg_malloc(framesize*sizeof(float));
 253   pcm[0]=_ogg_malloc(framesize*sizeof(float));
 254   pcm[1]=_ogg_malloc(framesize*sizeof(float));
 255   out[0]=_ogg_calloc(framesize/2,sizeof(float));
 256   out[1]=_ogg_calloc(framesize/2,sizeof(float));
 257   flr=_ogg_malloc(framesize*sizeof(float));
 258   lpc=_ogg_malloc(order*sizeof(float));
 259   buffer=_ogg_malloc(framesize*4);
 260   buffer2=buffer+framesize*2;
 261   window=_vorbis_window(0,framesize,framesize/2,framesize/2);
 262   mdct_init(&m_look,framesize);
 263   drft_init(&f_look,framesize);
 264   _vp_psy_init(&p_look,&_psy_set0,framesize/2,44100);
 265   floorinit(&floorlook,framesize/2,order,map);
 266
 267   for(i=0;i<P_BANDS;i++)
 268     for(j=0;j<P_LEVELS;j++)
 269       analysis("Ptonecurve",i*100+j,p_look.tonecurves[i][j],EHMER_MAX,0,0);
 270
 271   /* we cheat on the WAV header; we just bypass 44 bytes and never
 272      verify that it matches 16bit/stereo/44.1kHz. */
 273
 274   fread(buffer,1,44,stdin);
 275   fwrite(buffer,1,44,stdout);
 276   memset(buffer,0,framesize*2);
 277
 278   analysis("window",0,window,framesize,0,0);
 279
 280   fprintf(stderr,"Processing for frame size %d...\n",framesize);
 281
 282   while(!eos){
 283     long bytes=fread(buffer2,1,framesize*2,stdin);
 284     if(bytes<framesize*2)
 285       memset(buffer2+bytes,0,framesize*2-bytes);
 286
 287     if(bytes!=0){
 288
 289       /* uninterleave samples */
 290       for(i=0;i<framesize;i++){
 291         pcm[0][i]=((buffer[i*4+1]<<8)|
 292                       (0x00ff&(int)buffer[i*4]))/32768.f;
 293         pcm[1][i]=((buffer[i*4+3]<<8)|
 294                    (0x00ff&(int)buffer[i*4+2]))/32768.f;
 295       }
 296
 297       {
 298         float secs=framesize/44100.;
 299
 300         ampmax+=secs*ampmax_att_per_sec;
 301         if(ampmax<-9999)ampmax=-9999;
 302       }
 303       newmax=ampmax;
 304
 305       for(i=0;i<2;i++){
 306         float amp;
 307
 308         analysis("pre",frameno,pcm[i],framesize,0,0);
 309         memcpy(mask,pcm[i],sizeof(float)*framesize);
 310
 311         /* do the psychacoustics */
 312         for(j=0;j<framesize;j++)
 313           mask[j]=pcm[i][j]*=window[j];
 314
 315         drft_forward(&f_look,mask);
 316
 317         mask[0]/=(framesize/4.);
 318         for(j=1;j<framesize-1;j+=2)
 319           mask[(j+1)>>1]=4*hypot(mask[j],mask[j+1])/framesize;
 320
 321         mdct_forward(&m_look,pcm[i],pcm[i]);
 322         memcpy(mask+framesize/2,pcm[i],sizeof(float)*framesize/2);
 323         analysis("mdct",frameno,pcm[i],framesize/2,0,1);
 324         analysis("fft",frameno,mask,framesize/2,0,1);
 325
 326         {
 327           float ret;
 328           ret=_vp_compute_mask(&p_look,mask,mask+framesize/2,flr,NULL,ampmax);
 329           if(ret>newmax)newmax=ret;
 330         }
 331
 332         analysis("mask",frameno,flr,framesize/2,0,0);
 333
 334         for(j=0;j<framesize/2;j++)
 335           mask[j]=flr[j]+140.;
 336
 337         amp=sqrt(_curve_to_lpc(mask,mask,&floorlook));
 338         vorbis_lpc_to_lsp(mask,mask,floorlook.m);
 339         vorbis_lsp_to_curve(flr,floorlook.linearmap,floorlook.n,floorlook.ln,
 340                             mask,floorlook.m,amp,140.);
 341
 342         analysis("floor",frameno,flr,framesize/2,0,1);
 343
 344         _vp_apply_floor(&p_look,pcm[i],flr);
 345
 346
 347         analysis("quant",frameno,pcm[i],framesize/2,0,0);
 348
 349         /* re-add floor */
 350         for(j=0;j<framesize/2;j++){
 351           float val=rint(pcm[i][j]);
 352           tot++;
 353           if(val){
 354             nonz++;
 355             acc+=log(fabs(val)*2.f+1.f)/log(2);
 356             pcm[i][j]=val*flr[j];
 357           }else{
 358             pcm[i][j]=0.f;
 359           }
 360         }
 361
 362         analysis("final",frameno,pcm[i],framesize/2,0,1);
 363
 364         /* take it back to time */
 365         mdct_backward(&m_look,pcm[i],pcm[i]);
 366         for(j=0;j<framesize/2;j++)
 367           out[i][j]+=pcm[i][j]*window[j];
 368
 369         frameno++;
 370       }
 371       ampmax=newmax;
 372
 373       /* write data.  Use the part of buffer we're about to shift out */
 374       for(i=0;i<2;i++){
 375         char  *ptr=buffer+i*2;
 376         float *mono=out[i];
 377         for(j=0;j<framesize/2;j++){
 378           int val=mono[j]*32767.;
 379           /* might as well guard against clipping */
 380           if(val>32767)val=32767;
 381           if(val<-32768)val=-32768;
 382           ptr[0]=val&0xff;
 383           ptr[1]=(val>>8)&0xff;
 384           ptr+=4;
 385         }
 386       }
 387
 388       fprintf(stderr,"*");
 389       fwrite(buffer,1,framesize*2,stdout);
 390       memmove(buffer,buffer2,framesize*2);
 391
 392       for(i=0;i<2;i++){
 393         for(j=0,k=framesize/2;j<framesize/2;j++,k++)
 394           out[i][j]=pcm[i][k]*window[k];
 395       }
 396     }else
 397       eos=1;
 398   }
 399   fprintf(stderr,"average raw bits of entropy: %.03g/sample\n",acc/tot);
 400   fprintf(stderr,"average nonzero samples: %.03g/%d\n",nonz/tot*framesize/2,
 401           framesize/2);
 402   fprintf(stderr,"Done\n\n");
 403   return 0;
 404 }