lib/psytune.c

   1 /********************************************************************
   2  *                                                                  *
   3  * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE.   *
   4  * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
   5  * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
   6  * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
   7  *                                                                  *
   8  * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001             *
   9  * by the XIPHOPHORUS Company http://www.xiph.org/                  *
  10
  11  ********************************************************************
  12
  13  function: simple utility that runs audio through the psychoacoustics
  14            without encoding
  15  last mod: $Id: psytune.c,v 1.15 2001/05/27 06:44:00 xiphmont Exp $
  16
  17  ********************************************************************/
  18
  19 #include <stdio.h>
  20 #include <stdlib.h>
  21 #include <string.h>
  22 #include <math.h>
  23
  24 #include "vorbis/codec.h"
  25 #include "os.h"
  26 #include "psy.h"
  27 #include "mdct.h"
  28 #include "smallft.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
  38   -100.f,
  39   -140.f,
  40
  41   8,
  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    1,/* tonemaskp */
  47   /*  0   10   20   30   40   50   60   70   80   90   100 */
  48   {
  49    {-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f}, /*63*/
  50    {-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f}, /*88*/
  51    {-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f,-999.f}, /*125*/
  52
  53    {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*175*/
  54    {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*250*/
  55    {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*350*/
  56    {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*500*/
  57    {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*700*/
  58
  59    {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*1000*/
  60    {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*1400*/
  61    {-40.f,-40.f,-40.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*2000*/
  62    {-40.f,-40.f,-40.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*2800*/
  63    {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*4000*/
  64    {-35.f,-35.f,-35.f,-40.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*5600*/
  65
  66    {-30.f,-30.f,-33.f,-35.f,-40.f,-50.f,-60.f,-70.f,-80.f,-90.f,-100.f}, /*8000*/
  67    {-30.f,-30.f,-33.f,-35.f,-35.f,-45.f,-50.f,-60.f,-70.f,-90.f,-100.f}, /*11500*/
  68    {-24.f,-24.f,-26.f,-32.f,-32.f,-42.f,-50.f,-60.f,-70.f,-90.f,-100.f}, /*16000*/
  69
  70   },
  71
  72   1,/* peakattp */
  73   {{-14.f,-16.f,-18.f,-19.f,-20.f,-21.f,-22.f,-22.f,-24.f,-24.f,-24.f},/*63*/
  74    {-14.f,-16.f,-18.f,-19.f,-20.f,-21.f,-22.f,-22.f,-24.f,-24.f,-24.f},/*88*/
  75    {-14.f,-16.f,-18.f,-19.f,-20.f,-21.f,-22.f,-22.f,-24.f,-24.f,-24.f},/*125*/
  76    {-10.f,-10.f,-10.f,-10.f,-16.f,-16.f,-18.f,-20.f,-24.f,-24.f,-24.f},/*175*/
  77    {-10.f,-10.f,-10.f,-10.f,-16.f,-16.f,-18.f,-20.f,-24.f,-24.f,-24.f},/*250*/
  78    {-10.f,-10.f,-10.f,-10.f,-16.f,-16.f,-18.f,-20.f,-22.f,-24.f,-24.f},/*350*/
  79    {-10.f,-10.f,-10.f,-10.f,-16.f,-16.f,-18.f,-20.f,-22.f,-24.f,-24.f},/*500*/
  80    {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*700*/
  81    {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*1000*/
  82    {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*1400*/
  83    {-10.f,-10.f,-10.f,-10.f,-14.f,-14.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*2000*/
  84    {-10.f,-10.f,-10.f,-12.f,-16.f,-16.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*2400*/
  85    {-10.f,-10.f,-10.f,-12.f,-16.f,-16.f,-16.f,-20.f,-22.f,-24.f,-24.f},/*4000*/
  86    {-10.f,-10.f,-10.f,-12.f,-12.f,-14.f,-16.f,-18.f,-22.f,-24.f,-24.f},/*5600*/
  87    {-10.f,-10.f,-10.f,-10.f,-10.f,-14.f,-16.f,-18.f,-22.f,-24.f,-24.f},/*8000*/
  88    {-10.f,-10.f,-10.f,-10.f,-10.f,-14.f,-16.f,-18.f,-22.f,-24.f,-24.f},/*11500*/
  89    {-10.f,-10.f,-10.f,-10.f,-10.f,-12.f,-16.f,-18.f,-22.f,-24.f,-24.f},/*16000*/
  90   },
  91
  92   1,/*noisemaskp */
  93   -24.f,  /* suppress any noise curve over maxspec+n */
  94   .5f,   /* low window */
  95   .5f,   /* high window */
  96   25,
  97   25,
  98   {.000f, 0.f, /*63*/
  99    .000f, 0.f, /*88*/
 100    .000f, 0.f, /*125*/
 101    .000f, 0.f, /*175*/
 102    .000f, 0.f, /*250*/
 103    .000f, 0.f, /*350*/
 104    .000f, 0.f, /*500*/
 105    .200f, 0.f, /*700*/
 106    .300f, 0.f, /*1000*/
 107    .400f, 0.f, /*1400*/
 108    .400f, 0.f, /*2000*/
 109    .400f, 0.f, /*2800*/
 110    .700f, 0.f, /*4000*/
 111    .850f, 0.f, /*5600*/
 112    .900f, 0.f, /*8000*/
 113    .900f, 0.f, /*11500*/
 114    .900f, 1.f, /*16000*/
 115   },
 116
 117   95.f,  /* even decade + 5 is important; saves an rint() later in a
 118             tight loop) */
 119   -28.,
 120
 121 };
 122
 123 static int noisy=1;
 124 void analysis(char *base,int i,float *v,int n,int bark,int dB){
 125   if(noisy){
 126     int j;
 127     FILE *of;
 128     char buffer[80];
 129     sprintf(buffer,"%s_%d.m",base,i);
 130     of=fopen(buffer,"w");
 131
 132     for(j=0;j<n;j++){
 133       if(dB && v[j]==0)
 134           fprintf(of,"\n\n");
 135       else{
 136         if(bark)
 137           fprintf(of,"%g ",toBARK(22050.f*j/n));
 138         else
 139           fprintf(of,"%g ",(float)j);
 140
 141         if(dB){
 142           fprintf(of,"%g\n",todB(fabs(v+j)));
 143         }else{
 144           fprintf(of,"%g\n",v[j]);
 145         }
 146       }
 147     }
 148     fclose(of);
 149   }
 150 }
 151
 152 long frameno=0;
 153
 154 /****************************************************************/
 155
 156 int main(int argc,char *argv[]){
 157   int eos=0;
 158   float nonz=0.f;
 159   float acc=0.f;
 160   float tot=0.f;
 161   float ampmax=-9999,newmax;
 162
 163   int framesize=2048;
 164   int order=30;
 165   int map=256;
 166   float ampmax_att_per_sec=-10.;
 167
 168   float *pcm[2],*out[2],*window,*lpc,*flr,*mask;
 169   signed char *buffer,*buffer2;
 170   mdct_lookup m_look;
 171   drft_lookup f_look;
 172   drft_lookup f_look2;
 173   vorbis_look_psy p_look;
 174   long i,j,k;
 175
 176   int ath=0;
 177   int decayp=0;
 178
 179   argv++;
 180   while(*argv){
 181     if(*argv[0]=='-'){
 182       /* option */
 183       if(argv[0][1]=='v'){
 184         noisy=0;
 185       }
 186       if(argv[0][1]=='o'){
 187         order=atoi(argv[0]+2);
 188       }
 189       if(argv[0][1]=='m'){
 190         map=atoi(argv[0]+2);
 191       }
 192     }else
 193       if(*argv[0]=='+'){
 194         /* option */
 195         if(argv[0][1]=='v'){
 196           noisy=1;
 197         }
 198         if(argv[0][1]=='o'){
 199           order=atoi(argv[0]+2);
 200         }
 201         if(argv[0][1]=='m'){
 202           map=atoi(argv[0]+2);
 203         }
 204       }else
 205         framesize=atoi(argv[0]);
 206     argv++;
 207   }
 208
 209   mask=_ogg_malloc(framesize*sizeof(float));
 210   pcm[0]=_ogg_malloc(framesize*sizeof(float));
 211   pcm[1]=_ogg_malloc(framesize*sizeof(float));
 212   out[0]=_ogg_calloc(framesize/2,sizeof(float));
 213   out[1]=_ogg_calloc(framesize/2,sizeof(float));
 214   flr=_ogg_malloc(framesize*sizeof(float));
 215   lpc=_ogg_malloc(order*sizeof(float));
 216   buffer=_ogg_malloc(framesize*4);
 217   buffer2=buffer+framesize*2;
 218   window=_vorbis_window(0,framesize,framesize/2,framesize/2);
 219   mdct_init(&m_look,framesize);
 220   drft_init(&f_look,framesize);
 221   drft_init(&f_look2,framesize/2);
 222   _vp_psy_init(&p_look,&_psy_set0,framesize/2,44100);
 223
 224   for(i=0;i<P_BANDS;i++)
 225     for(j=0;j<P_LEVELS;j++)
 226       analysis("Ptonecurve",i*100+j,p_look.tonecurves[i][j],EHMER_MAX,0,0);
 227
 228   /* we cheat on the WAV header; we just bypass 44 bytes and never
 229      verify that it matches 16bit/stereo/44.1kHz. */
 230
 231   fread(buffer,1,44,stdin);
 232   fwrite(buffer,1,44,stdout);
 233   memset(buffer,0,framesize*2);
 234
 235   analysis("window",0,window,framesize,0,0);
 236
 237   fprintf(stderr,"Processing for frame size %d...\n",framesize);
 238
 239   while(!eos){
 240     long bytes=fread(buffer2,1,framesize*2,stdin);
 241     if(bytes<framesize*2)
 242       memset(buffer2+bytes,0,framesize*2-bytes);
 243
 244     if(bytes!=0){
 245
 246       /* uninterleave samples */
 247       for(i=0;i<framesize;i++){
 248         pcm[0][i]=((buffer[i*4+1]<<8)|
 249                       (0x00ff&(int)buffer[i*4]))/32768.f;
 250         pcm[1][i]=((buffer[i*4+3]<<8)|
 251                    (0x00ff&(int)buffer[i*4+2]))/32768.f;
 252       }
 253
 254       {
 255         float secs=framesize/44100.;
 256
 257         ampmax+=secs*ampmax_att_per_sec;
 258         if(ampmax<-9999)ampmax=-9999;
 259       }
 260       newmax=ampmax;
 261
 262       for(i=0;i<2;i++){
 263         float amp;
 264
 265         analysis("pre",frameno,pcm[i],framesize,0,0);
 266         memcpy(mask,pcm[i],sizeof(float)*framesize);
 267
 268         /* do the psychacoustics */
 269         for(j=0;j<framesize;j++)
 270           mask[j]=pcm[i][j]*=window[j];
 271
 272         drft_forward(&f_look,mask);
 273
 274         mask[0]/=(framesize/4.);
 275         for(j=1;j<framesize-1;j+=2)
 276           mask[(j+1)>>1]=4*hypot(mask[j],mask[j+1])/framesize;
 277
 278         mdct_forward(&m_look,pcm[i],pcm[i]);
 279         memcpy(mask+framesize/2,pcm[i],sizeof(float)*framesize/2);
 280         analysis("mdct",frameno,pcm[i],framesize/2,0,1);
 281         analysis("fft",frameno,mask,framesize/2,0,1);
 282
 283         {
 284           float ret;
 285           ret=_vp_compute_mask(&p_look,mask,mask+framesize/2,flr,NULL,ampmax);
 286           if(ret>newmax)newmax=ret;
 287         }
 288
 289         analysis("mask",frameno,flr,framesize/2,0,0);
 290
 291         mask[framesize-1]=0.;
 292         mask[0]=0.;
 293         for(j=1;j<framesize-1;j+=2){
 294           mask[j]=todB(pcm[i]+((j+1)>>1));
 295           mask[j+1]=0;
 296         }
 297
 298         analysis("lfft",frameno,mask,framesize,0,0);
 299         drft_backward(&f_look,mask);
 300
 301         analysis("cep",frameno,mask,framesize,0,0);
 302         analysis("logcep",frameno,mask,framesize,0,1);
 303
 304
 305
 306         /*for(j=0;j<framesize/2;j++){
 307           float val=fromdB(flr[j]);
 308           int p=rint(pcm[i][j]/val);
 309           pcm[i][j]=p*val;
 310           }*/
 311
 312         /*for(j=0;j<framesize/2;j++){
 313           float val=todB(pcm[i]+j);
 314           if(val+6.<flr[j])
 315             pcm[i][j]=0.;
 316             }*/
 317
 318         for(j=0;j<framesize/2;j++){
 319           float val=rint(todB(pcm[i]+j)/6);
 320           if(pcm[i][j]>0)
 321             pcm[i][j]=fromdB(val*6);
 322           else
 323             pcm[i][j]=-fromdB(val*6);
 324         }
 325
 326
 327         analysis("final",frameno,pcm[i],framesize/2,0,1);
 328
 329         /* take it back to time */
 330         mdct_backward(&m_look,pcm[i],pcm[i]);
 331         for(j=0;j<framesize/2;j++)
 332           out[i][j]+=pcm[i][j]*window[j];
 333
 334         frameno++;
 335       }
 336       ampmax=newmax;
 337
 338       /* write data.  Use the part of buffer we're about to shift out */
 339       for(i=0;i<2;i++){
 340         char  *ptr=buffer+i*2;
 341         float *mono=out[i];
 342         for(j=0;j<framesize/2;j++){
 343           int val=mono[j]*32767.;
 344           /* might as well guard against clipping */
 345           if(val>32767)val=32767;
 346           if(val<-32768)val=-32768;
 347           ptr[0]=val&0xff;
 348           ptr[1]=(val>>8)&0xff;
 349           ptr+=4;
 350         }
 351       }
 352
 353       fprintf(stderr,"*");
 354       fwrite(buffer,1,framesize*2,stdout);
 355       memmove(buffer,buffer2,framesize*2);
 356
 357       for(i=0;i<2;i++){
 358         for(j=0,k=framesize/2;j<framesize/2;j++,k++)
 359           out[i][j]=pcm[i][k]*window[k];
 360       }
 361     }else
 362       eos=1;
 363   }
 364   fprintf(stderr,"average raw bits of entropy: %.03g/sample\n",acc/tot);
 365   fprintf(stderr,"average nonzero samples: %.03g/%d\n",nonz/tot*framesize/2,
 366           framesize/2);
 367   fprintf(stderr,"Done\n\n");
 368   return 0;
 369 }