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