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