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