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