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.8 2000/11/06 00:07:02 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 "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   1,/*decayp*/
  37   1,/*smoothp*/
  38   0,.2,
  39
  40   -100.,
  41   -140.,
  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   {{-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*63*/
  50    {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*88*/
  51    {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*125*/
  52    {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*175*/
  53    {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*250*/
  54    {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*350*/
  55    {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*500*/
  56    {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*700*/
  57    {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*1000*/
  58    {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*1400*/
  59    {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*2000*/
  60    {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*2800*/
  61    {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*4000*/
  62    {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*5600*/
  63    {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*8000*/
  64    {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*11500*/
  65    {-40.,-40.,-40.,-40.,-40.,-50.,-60.,-70.,-80.,-90.,-100.}, /*16000*/
  66   },
  67
  68   1,/* peakattp */
  69   {{-14.,-16.,-18.,-19.,-20.,-21.,-22.,-22.,-22.,-22.,-22.}, /*63*/
  70    {-14.,-16.,-18.,-19.,-20.,-21.,-22.,-22.,-22.,-22.,-22.}, /*88*/
  71    {-14.,-16.,-18.,-19.,-20.,-21.,-22.,-22.,-22.,-22.,-22.}, /*125*/
  72    { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  73    { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  74    { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  75    { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  76    { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  77    { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  78    { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  79    { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  80    { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  81    { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  82    { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-18.,-20.}, /*8000*/
  83    { -7., -8., -9.,-10.,-10.,-11.,-12.,-13.,-15.,-16.,-17.}, /*8000*/
  84    { -6., -7., -9., -9., -9., -9.,-10.,-11.,-12.,-13.,-14.}, /*11500*/
  85    { -6., -6., -9., -9., -9., -9., -9., -9.,-10.,-11.,-12.}, /*16000*/
  86   },
  87
  88   1,/*noisemaskp */
  89   /*  0   10   20   30   40   50   60    70    80    90   100 */
  90   {{-99.,-99.,-99.,-99.,-99.,-99.,-99.,-199.,-199.,-199.,-199.}, /*63*/
  91    {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-199.,-199.,-199.,-199.}, /*88*/
  92    {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-199.,-199.,-199.,-199.}, /*125*/
  93    {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-199.,-199.,-199.,-199.}, /*175*/
  94    {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-199.,-199.,-199.,-199.}, /*250*/
  95    {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-199.,-199.,-199.,-199.}, /*350*/
  96    {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-199.,-199.,-199.,-199.}, /*350*/
  97    {-99.,-99.,-99.,-99.,-99.,-99.,-99.,-199.,-199.,-199.,-199.}, /*350*/
  98    { -6., -6., -6., -6., -6., -6., -6.,  -6.,  -6.,  -6.,  -6.}, /*2000*/
  99    { -6., -6., -6., -6., -6., -6., -6.,  -6.,  -6.,  -6.,  -6.}, /*2000*/
 100    { -6., -6., -6., -6., -6., -6., -6.,  -6.,  -6.,  -6.,  -6.}, /*2000*/
 101    { -6., -6., -6., -6., -6., -6., -6.,  -6.,  -6.,  -6.,  -6.}, /*2800*/
 102    { -0., -0., -0., -0., -0., -0., -0.,  -0.,  -0.,  -0.,  -0.}, /*4000*/
 103    {  0.,  0.,  0.,  0.,  0.,  0.,  0.,   0.,   0.,   0.,   0.}, /*5600*/
 104    {  4.,  4.,  4.,  4.,  4.,  4.,  4.,   4.,   4.,   4.,   4.}, /*8000*/
 105    {  4.,  4.,  4.,  4.,  4.,  4.,  4.,   4.,   4.,   4.,   4.}, /*11500*/
 106    {  4.,  4.,  4.,  4.,  4.,  4.,  4.,   4.,   4.,   4.,   4.}, /*16000*/
 107   },
 108
 109   110.,
 110
 111   -0, -.005,  /* attack/decay control */
 112 };
 113
 114 static int noisy=0;
 115 void analysis(char *base,int i,float *v,int n,int bark,int dB){
 116   if(noisy){
 117     int j;
 118     FILE *of;
 119     char buffer[80];
 120     sprintf(buffer,"%s_%d.m",base,i);
 121     of=fopen(buffer,"w");
 122
 123     for(j=0;j<n;j++){
 124       if(dB && v[j]==0)
 125           fprintf(of,"\n\n");
 126       else{
 127         if(bark)
 128           fprintf(of,"%g ",toBARK(22050.*j/n));
 129         else
 130           fprintf(of,"%g ",(float)j);
 131
 132         if(dB){
 133           fprintf(of,"%g\n",todB(fabs(v[j])));
 134         }else{
 135           fprintf(of,"%g\n",v[j]);
 136         }
 137       }
 138     }
 139     fclose(of);
 140   }
 141 }
 142
 143 typedef struct {
 144   long n;
 145   int ln;
 146   int  m;
 147   int *linearmap;
 148
 149   vorbis_info_floor0 *vi;
 150   lpc_lookup lpclook;
 151 } vorbis_look_floor0;
 152
 153 long granulepos=0;
 154
 155 /* hacked from floor0.c */
 156 static void floorinit(vorbis_look_floor0 *look,int n,int m,int ln){
 157   int j;
 158   float scale;
 159   look->m=m;
 160   look->n=n;
 161   look->ln=ln;
 162   lpc_init(&look->lpclook,look->ln,look->m);
 163
 164   scale=look->ln/toBARK(22050.);
 165
 166   look->linearmap=_ogg_malloc(look->n*sizeof(int));
 167   for(j=0;j<look->n;j++){
 168     int val=floor( toBARK(22050./n*j) *scale);
 169     if(val>look->ln)val=look->ln;
 170     look->linearmap[j]=val;
 171   }
 172 }
 173
 174 int main(int argc,char *argv[]){
 175   int eos=0;
 176   float nonz=0.;
 177   float acc=0.;
 178   float tot=0.;
 179
 180   int framesize=2048;
 181   int order=32;
 182   int map=256;
 183
 184   float *pcm[2],*out[2],*window,*decay[2],*lpc,*floor;
 185   signed char *buffer,*buffer2;
 186   mdct_lookup m_look;
 187   vorbis_look_psy p_look;
 188   long i,j,k;
 189
 190   vorbis_look_floor0 floorlook;
 191
 192   int ath=0;
 193   int decayp=0;
 194
 195   argv++;
 196   while(*argv){
 197     if(*argv[0]=='-'){
 198       /* option */
 199       if(argv[0][1]=='v'){
 200         noisy=0;
 201       }
 202       if(argv[0][1]=='A'){
 203         ath=0;
 204       }
 205       if(argv[0][1]=='D'){
 206         decayp=0;
 207       }
 208       if(argv[0][1]=='X'){
 209         ath=0;
 210         decayp=0;
 211       }
 212     }else
 213       if(*argv[0]=='+'){
 214         /* option */
 215         if(argv[0][1]=='v'){
 216           noisy=1;
 217         }
 218         if(argv[0][1]=='A'){
 219           ath=1;
 220         }
 221         if(argv[0][1]=='D'){
 222           decayp=1;
 223         }
 224         if(argv[0][1]=='X'){
 225           ath=1;
 226           decayp=1;
 227         }
 228       }else
 229         framesize=atoi(argv[0]);
 230     argv++;
 231   }
 232
 233   pcm[0]=_ogg_malloc(framesize*sizeof(float));
 234   pcm[1]=_ogg_malloc(framesize*sizeof(float));
 235   out[0]=_ogg_calloc(framesize/2,sizeof(float));
 236   out[1]=_ogg_calloc(framesize/2,sizeof(float));
 237   decay[0]=_ogg_calloc(framesize/2,sizeof(float));
 238   decay[1]=_ogg_calloc(framesize/2,sizeof(float));
 239   floor=_ogg_malloc(framesize*sizeof(float));
 240   lpc=_ogg_malloc(order*sizeof(float));
 241   buffer=_ogg_malloc(framesize*4);
 242   buffer2=buffer+framesize*2;
 243   window=_vorbis_window(0,framesize,framesize/2,framesize/2);
 244   mdct_init(&m_look,framesize);
 245   _vp_psy_init(&p_look,&_psy_set0,framesize/2,44100);
 246   floorinit(&floorlook,framesize/2,order,map);
 247
 248   for(i=0;i<P_BANDS;i++)
 249     for(j=0;j<P_LEVELS;j++)
 250       analysis("Ptonecurve",i*100+j,p_look.tonecurves[i][j],EHMER_MAX,0,1);
 251
 252   /* we cheat on the WAV header; we just bypass 44 bytes and never
 253      verify that it matches 16bit/stereo/44.1kHz. */
 254
 255   fread(buffer,1,44,stdin);
 256   fwrite(buffer,1,44,stdout);
 257   memset(buffer,0,framesize*2);
 258
 259   analysis("window",0,window,framesize,0,0);
 260
 261   fprintf(stderr,"Processing for frame size %d...\n",framesize);
 262
 263   while(!eos){
 264     long bytes=fread(buffer2,1,framesize*2,stdin);
 265     if(bytes<framesize*2)
 266       memset(buffer2+bytes,0,framesize*2-bytes);
 267
 268     if(bytes!=0){
 269
 270       /* uninterleave samples */
 271       for(i=0;i<framesize;i++){
 272         pcm[0][i]=((buffer[i*4+1]<<8)|
 273                       (0x00ff&(int)buffer[i*4]))/32768.;
 274         pcm[1][i]=((buffer[i*4+3]<<8)|
 275                    (0x00ff&(int)buffer[i*4+2]))/32768.;
 276       }
 277
 278       for(i=0;i<2;i++){
 279         float amp;
 280
 281         analysis("pre",granulepos,pcm[i],framesize,0,0);
 282
 283         /* do the psychacoustics */
 284         for(j=0;j<framesize;j++)
 285           pcm[i][j]*=window[j];
 286
 287         mdct_forward(&m_look,pcm[i],pcm[i]);
 288
 289         analysis("mdct",granulepos,pcm[i],framesize/2,1,1);
 290
 291         _vp_compute_mask(&p_look,pcm[i],floor,decay[i]);
 292
 293         analysis("floor",frameno,floor,framesize/2,1,1);
 294         analysis("decay",frameno,decay[i],framesize/2,1,1);
 295
 296         _vp_apply_floor(&p_look,pcm[i],floor);
 297         /*r(j=0;j<framesize/2;j++)
 298           if(fabs(pcm[i][j])<1.)pcm[i][j]=0;*/
 299
 300         analysis("quant",granulepos,pcm[i],framesize/2,1,1);
 301
 302         /* re-add floor */
 303         for(j=0;j<framesize/2;j++){
 304           float val=rint(pcm[i][j]);
 305           tot++;
 306           if(val){
 307             nonz++;
 308             acc+=log(fabs(val)*2.+1.)/log(2);
 309             pcm[i][j]=val*floor[j];
 310           }else{
 311             pcm[i][j]=0;
 312           }
 313         }
 314
 315         analysis("final",granulepos,pcm[i],framesize/2,1,1);
 316
 317         /* take it back to time */
 318         mdct_backward(&m_look,pcm[i],pcm[i]);
 319         for(j=0;j<framesize/2;j++)
 320           out[i][j]+=pcm[i][j]*window[j];
 321
 322         granulepos++;
 323       }
 324
 325       /* write data.  Use the part of buffer we're about to shift out */
 326       for(i=0;i<2;i++){
 327         char  *ptr=buffer+i*2;
 328         float *mono=out[i];
 329         for(j=0;j<framesize/2;j++){
 330           int val=mono[j]*32767.;
 331           /* might as well guard against clipping */
 332           if(val>32767)val=32767;
 333           if(val<-32768)val=-32768;
 334           ptr[0]=val&0xff;
 335           ptr[1]=(val>>8)&0xff;
 336           ptr+=4;
 337         }
 338       }
 339
 340       fprintf(stderr,"*");
 341       fwrite(buffer,1,framesize*2,stdout);
 342       memmove(buffer,buffer2,framesize*2);
 343
 344       for(i=0;i<2;i++){
 345         for(j=0,k=framesize/2;j<framesize/2;j++,k++)
 346           out[i][j]=pcm[i][k]*window[k];
 347       }
 348     }else
 349       eos=1;
 350   }
 351   fprintf(stderr,"average raw bits of entropy: %.03g/sample\n",acc/tot);
 352   fprintf(stderr,"average nonzero samples: %.03g/%d\n",nonz/tot*framesize/2,
 353           framesize/2);
 354   fprintf(stderr,"Done\n\n");
 355   return 0;
 356 }