lib/psytune.c

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