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