Cascading fully functional
[platform/upstream/libvorbis.git] / lib / res0.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: residue backend 0 and 1 implementation
14  last mod: $Id: res0.c,v 1.28 2001/06/04 05:50:10 xiphmont Exp $
15
16  ********************************************************************/
17
18 /* Slow, slow, slow, simpleminded and did I mention it was slow?  The
19    encode/decode loops are coded for clarity and performance is not
20    yet even a nagging little idea lurking in the shadows.  Oh and BTW,
21    it's slow. */
22
23 #include <stdlib.h>
24 #include <string.h>
25 #include <math.h>
26 #include <stdio.h>
27 #include <ogg/ogg.h>
28 #include "vorbis/codec.h"
29 #include "codec_internal.h"
30 #include "registry.h"
31 #include "codebook.h"
32 #include "misc.h"
33 #include "os.h"
34
35 typedef struct {
36   vorbis_info_residue0 *info;
37   int         map;
38   
39   int         parts;
40   int         stages;
41   codebook   *fullbooks;
42   codebook   *phrasebook;
43   codebook ***partbooks;
44
45   int         partvals;
46   int       **decodemap;
47
48   long        phrase;
49   long        bits[32];
50   long        vals[32];
51 } vorbis_look_residue0;
52
53 vorbis_info_residue *res0_copy_info(vorbis_info_residue *vr){
54   vorbis_info_residue0 *info=(vorbis_info_residue0 *)vr;
55   vorbis_info_residue0 *ret=_ogg_malloc(sizeof(vorbis_info_residue0));
56   memcpy(ret,info,sizeof(vorbis_info_residue0));
57   return(ret);
58 }
59
60 void res0_free_info(vorbis_info_residue *i){
61   if(i){
62     memset(i,0,sizeof(vorbis_info_residue0));
63     _ogg_free(i);
64   }
65 }
66
67 void res0_free_look(vorbis_look_residue *i){
68   int j;
69   if(i){
70     vorbis_look_residue0 *look=(vorbis_look_residue0 *)i;
71     long resbitsT=0;
72     long resvalsT=0;
73
74     for(j=0;j<look->parts;j++)resbitsT+=look->bits[j];
75     for(j=0;j<look->parts;j++)resvalsT+=look->vals[j];
76     fprintf(stderr,
77             "Encoded %ld res vectors in %ld phrasing and %ld res bits\n\t",
78             resvalsT,look->phrase,resbitsT);
79     for(j=1;j<look->parts;j++)
80       fprintf(stderr,"%ld(%ld):%ld  %g/sample\n",j,look->vals[j],look->bits[j],
81               (float)look->bits[j]/look->vals[j]);
82     fprintf(stderr,"\n");
83  
84
85     for(j=0;j<look->parts;j++)
86       if(look->partbooks[j])_ogg_free(look->partbooks[j]);
87     _ogg_free(look->partbooks);
88     for(j=0;j<look->partvals;j++)
89       _ogg_free(look->decodemap[j]);
90     _ogg_free(look->decodemap);
91     memset(i,0,sizeof(vorbis_look_residue0));
92     _ogg_free(i);
93   }
94 }
95
96 static int ilog(unsigned int v){
97   int ret=0;
98   while(v){
99     ret++;
100     v>>=1;
101   }
102   return(ret);
103 }
104
105 static int icount(unsigned int v){
106   int ret=0;
107   while(v){
108     ret+=v&1;
109     v>>=1;
110   }
111   return(ret);
112 }
113
114
115 void res0_pack(vorbis_info_residue *vr,oggpack_buffer *opb){
116   vorbis_info_residue0 *info=(vorbis_info_residue0 *)vr;
117   int j,acc=0;
118   oggpack_write(opb,info->begin,24);
119   oggpack_write(opb,info->end,24);
120
121   oggpack_write(opb,info->grouping-1,24);  /* residue vectors to group and 
122                                              code with a partitioned book */
123   oggpack_write(opb,info->partitions-1,6); /* possible partition choices */
124   oggpack_write(opb,info->groupbook,8);  /* group huffman book */
125
126   /* secondstages is a bitmask; as encoding progresses pass by pass, a
127      bitmask of one indicates this partition class has bits to write
128      this pass */
129   for(j=0;j<info->partitions;j++){
130     if(ilog(info->secondstages[j])>3){
131       /* yes, this is a minor hack due to not thinking ahead */
132       oggpack_write(opb,info->secondstages[j],3); 
133       oggpack_write(opb,1,1);
134       oggpack_write(opb,info->secondstages[j]>>3,5); 
135     }else
136       oggpack_write(opb,info->secondstages[j],4); /* trailing zero */
137     acc+=icount(info->secondstages[j]);
138   }
139   for(j=0;j<acc;j++)
140     oggpack_write(opb,info->booklist[j],8);
141
142 }
143
144 /* vorbis_info is for range checking */
145 vorbis_info_residue *res0_unpack(vorbis_info *vi,oggpack_buffer *opb){
146   int j,acc=0;
147   vorbis_info_residue0 *info=_ogg_calloc(1,sizeof(vorbis_info_residue0));
148   codec_setup_info     *ci=vi->codec_setup;
149
150   info->begin=oggpack_read(opb,24);
151   info->end=oggpack_read(opb,24);
152   info->grouping=oggpack_read(opb,24)+1;
153   info->partitions=oggpack_read(opb,6)+1;
154   info->groupbook=oggpack_read(opb,8);
155
156   for(j=0;j<info->partitions;j++){
157     int cascade=info->secondstages[j]=oggpack_read(opb,3);
158     if(oggpack_read(opb,1))
159       cascade|=(oggpack_read(opb,5)<<3);
160     acc+=icount(cascade);
161   }
162   for(j=0;j<acc;j++)
163     info->booklist[j]=oggpack_read(opb,8);
164
165   if(info->groupbook>=ci->books)goto errout;
166   for(j=0;j<acc;j++)
167     if(info->booklist[j]>=ci->books)goto errout;
168
169   return(info);
170  errout:
171   res0_free_info(info);
172   return(NULL);
173 }
174
175 vorbis_look_residue *res0_look (vorbis_dsp_state *vd,vorbis_info_mode *vm,
176                           vorbis_info_residue *vr){
177   vorbis_info_residue0 *info=(vorbis_info_residue0 *)vr;
178   vorbis_look_residue0 *look=_ogg_calloc(1,sizeof(vorbis_look_residue0));
179   backend_lookup_state *be=vd->backend_state;
180
181   int j,k,acc=0;
182   int dim;
183   int maxstage=0;
184   look->info=info;
185   look->map=vm->mapping;
186
187   look->parts=info->partitions;
188   look->fullbooks=be->fullbooks;
189   look->phrasebook=be->fullbooks+info->groupbook;
190   dim=look->phrasebook->dim;
191
192   look->partbooks=_ogg_calloc(look->parts,sizeof(codebook **));
193
194   for(j=0;j<look->parts;j++){
195     int stages=ilog(info->secondstages[j]);
196     if(stages){
197       if(stages>maxstage)maxstage=stages;
198       look->partbooks[j]=_ogg_calloc(stages,sizeof(codebook *));
199       for(k=0;k<stages;k++)
200         if(info->secondstages[j]&(1<<k))
201           look->partbooks[j][k]=be->fullbooks+info->booklist[acc++];
202     }
203   }
204
205   look->partvals=rint(pow(look->parts,dim));
206   look->stages=maxstage;
207   look->decodemap=_ogg_malloc(look->partvals*sizeof(int *));
208   for(j=0;j<look->partvals;j++){
209     long val=j;
210     long mult=look->partvals/look->parts;
211     look->decodemap[j]=_ogg_malloc(dim*sizeof(int));
212     for(k=0;k<dim;k++){
213       long deco=val/mult;
214       val-=deco*mult;
215       mult/=look->parts;
216       look->decodemap[j][k]=deco;
217     }
218   }
219
220   return(look);
221 }
222
223
224 /* does not guard against invalid settings; eg, a subn of 16 and a
225    subgroup request of 32.  Max subn of 128 */
226 static int _interleaved_testhack(float *vec,int n,vorbis_look_residue0 *look,
227                                  int auxparts,int auxpartnum){
228   vorbis_info_residue0 *info=look->info;
229   int i,j=0;
230   float max,localmax=0.f;
231   float temp[128];
232   float entropy[8];
233
234   /* setup */
235   for(i=0;i<n;i++)temp[i]=fabs(vec[i]);
236
237   /* handle case subgrp==1 outside */
238   for(i=0;i<n;i++)
239     if(temp[i]>localmax)localmax=temp[i];
240   max=localmax;
241
242   for(i=0;i<n;i++)temp[i]=rint(temp[i]);
243   
244   while(1){
245     entropy[j]=localmax;
246     n>>=1;
247     if(!n)break;
248     j++;
249
250     for(i=0;i<n;i++){
251       temp[i]+=temp[i+n];
252     }
253     localmax=0.f;
254     for(i=0;i<n;i++)
255       if(temp[i]>localmax)localmax=temp[i];
256   }
257
258   for(i=0;i<auxparts-1;i++)
259     if(auxpartnum<info->blimit[i] &&
260        entropy[info->subgrp[i]]<=info->entmax[i] &&
261        max<=info->ampmax[i])
262       break;
263
264   return(i);
265 }
266
267 static int _testhack(float *vec,int n,vorbis_look_residue0 *look,
268                      int auxparts,int auxpartnum){
269   vorbis_info_residue0 *info=look->info;
270   int i,j=0;
271   float max,localmax=0.f;
272   float temp[128];
273   float entropy[8];
274
275   /* setup */
276   for(i=0;i<n;i++)temp[i]=fabs(vec[i]);
277
278   /* handle case subgrp==1 outside */
279   for(i=0;i<n;i++)
280     if(temp[i]>localmax)localmax=temp[i];
281   max=localmax;
282
283   for(i=0;i<n;i++)temp[i]=rint(temp[i]);
284   
285   while(n){
286     entropy[j]=localmax;
287     n>>=1;
288     j++;
289     if(!n)break;
290     for(i=0;i<n;i++){
291       temp[i]=temp[i*2]+temp[i*2+1];
292     }
293     localmax=0.f;
294     for(i=0;i<n;i++)
295       if(temp[i]>localmax)localmax=temp[i];
296   }
297
298   for(i=0;i<auxparts-1;i++)
299     if(auxpartnum<info->blimit[i] &&
300        entropy[info->subgrp[i]]<=info->entmax[i] &&
301        max<=info->ampmax[i])
302       break;
303
304   return(i);
305 }
306
307 static int _interleaved_encodepart(oggpack_buffer *opb,float *vec, int n,
308                                    codebook *book,vorbis_look_residue0 *look){
309   int i,bits=0;
310   int dim=book->dim;
311   int step=n/dim;
312 #ifdef TRAIN_RESENT      
313   char buf[80];
314   FILE *f;
315   sprintf(buf,"res0_b%d.vqd",book-look->fullbooks);
316   f=fopen(buf,"a");
317 #endif
318
319   for(i=0;i<step;i++){
320     int entry=vorbis_book_besterror(book,vec+i,step,0);
321
322 #ifdef TRAIN_RESENT      
323     fprintf(f,"%d\n",entry);
324 #endif
325
326     bits+=vorbis_book_encode(book,entry,opb);
327   }
328
329 #ifdef TRAIN_RESENT      
330   fclose(f);
331 #endif
332   return(bits);
333 }
334  
335 static int _encodepart(oggpack_buffer *opb,float *vec, int n,
336                        codebook *book,vorbis_look_residue0 *look){
337   int i,bits=0;
338   int dim=book->dim;
339   int step=n/dim;
340 #ifdef TRAIN_RESENT      
341   char buf[80];
342   FILE *f;
343   sprintf(buf,"res0_b%d.vqd",book-look->fullbooks);
344   f=fopen(buf,"a");
345 #endif
346
347   for(i=0;i<step;i++){
348     int entry=vorbis_book_besterror(book,vec+i*dim,1,0);
349
350 #ifdef TRAIN_RESENT      
351     fprintf(f,"%d\n",entry);
352 #endif
353
354     bits+=vorbis_book_encode(book,entry,opb);
355   }
356
357 #ifdef TRAIN_RESENT      
358   fclose(f);
359 #endif
360   return(bits);
361 }
362
363 static int _01forward(vorbis_block *vb,vorbis_look_residue *vl,
364                       float **in,int ch,
365                       int (*classify)(float *,int,vorbis_look_residue0 *,
366                                       int,int),
367                       int (*encode)(oggpack_buffer *,float *,int,
368                                     codebook *,vorbis_look_residue0 *)){
369   long i,j,k,l,s;
370   vorbis_look_residue0 *look=(vorbis_look_residue0 *)vl;
371   vorbis_info_residue0 *info=look->info;
372
373   /* move all this setup out later */
374   int samples_per_partition=info->grouping;
375   int possible_partitions=info->partitions;
376   int partitions_per_word=look->phrasebook->dim;
377   int n=info->end-info->begin;
378
379   int partvals=n/samples_per_partition;
380   int partwords=(partvals+partitions_per_word-1)/partitions_per_word;
381   long **partword=_vorbis_block_alloc(vb,ch*sizeof(long *));
382
383   partvals=partwords*partitions_per_word;
384
385   /* we find the patition type for each partition of each
386      channel.  We'll go back and do the interleaved encoding in a
387      bit.  For now, clarity */
388  
389   for(i=0;i<ch;i++){
390     partword[i]=_vorbis_block_alloc(vb,n/samples_per_partition*sizeof(long));
391     memset(partword[i],0,n/samples_per_partition*sizeof(long));
392   }
393
394   for(i=info->begin,l=0;i<info->end;i+=samples_per_partition,l++){
395     for(j=0;j<ch;j++)
396       /* do the partition decision based on the 'entropy'
397          int the block */
398       partword[j][l]=
399         classify(in[j]+i,samples_per_partition,look,possible_partitions,l);
400   
401   }
402
403   /* we code the partition words for each channel, then the residual
404      words for a partition per channel until we've written all the
405      residual words for that partition word.  Then write the next
406      partition channel words... */
407
408   for(s=0;s<look->stages;s++){
409     for(i=info->begin,l=0;i<info->end;){
410
411       /* first we encode a partition codeword for each channel */
412       if(s==0){
413         for(j=0;j<ch;j++){
414           long val=partword[j][l];
415           for(k=1;k<partitions_per_word;k++)
416             val= val*possible_partitions+partword[j][l+k];
417           look->phrase+=vorbis_book_encode(look->phrasebook,val,&vb->opb);
418         }
419       }
420
421       /* now we encode interleaved residual values for the partitions */
422       for(k=0;k<partitions_per_word;k++,l++,i+=samples_per_partition){
423         
424         for(j=0;j<ch;j++){
425           if(info->secondstages[partword[j][l]]&(1<<s)){
426             codebook *statebook=look->partbooks[partword[j][l]][s];
427             if(statebook){
428               int ret=encode(&vb->opb,in[j]+i,samples_per_partition,
429                                statebook,look);
430               look->bits[partword[j][l]]+=ret;
431               if(s==0)look->vals[partword[j][l]]+=samples_per_partition;
432             }
433           }
434         }
435       }
436     }
437   }
438
439   return(0);
440 }
441
442 /* a truncated packet here just means 'stop working'; it's not an error */
443 static int _01inverse(vorbis_block *vb,vorbis_look_residue *vl,
444                       float **in,int ch,
445                       long (*decodepart)(codebook *, float *, 
446                                          oggpack_buffer *,int,int)){
447
448   long i,j,k,l,s,transend=vb->pcmend/2;
449   vorbis_look_residue0 *look=(vorbis_look_residue0 *)vl;
450   vorbis_info_residue0 *info=look->info;
451
452   /* move all this setup out later */
453   int samples_per_partition=info->grouping;
454   int partitions_per_word=look->phrasebook->dim;
455   int n=info->end-info->begin;
456
457   int partvals=n/samples_per_partition;
458   int partwords=(partvals+partitions_per_word-1)/partitions_per_word;
459   int ***partword=alloca(ch*sizeof(int **));
460   float **work=alloca(ch*sizeof(float *));
461   partvals=partwords*partitions_per_word;
462
463   /* make sure we're zeroed up to the start */
464   for(j=0;j<ch;j++){
465     work[j]=_vorbis_block_alloc(vb,n*sizeof(float));
466     partword[j]=_vorbis_block_alloc(vb,partwords*sizeof(int *));
467     memset(work[j],0,sizeof(float)*n);
468   }
469
470   for(s=0;s<look->stages;s++){
471     for(i=info->begin,l=0;i<info->end;l++){
472
473       if(s==0){
474         /* fetch the partition word for each channel */
475         for(j=0;j<ch;j++){
476           int temp=vorbis_book_decode(look->phrasebook,&vb->opb);
477           if(temp==-1)goto eopbreak;
478           partword[j][l]=look->decodemap[temp];
479           if(partword[j][l]==NULL)goto errout;
480         }
481       }
482       
483       /* now we decode residual values for the partitions */
484       for(k=0;k<partitions_per_word;k++,i+=samples_per_partition)
485         for(j=0;j<ch;j++){
486           if(info->secondstages[partword[j][l][k]]&(1<<s)){
487             codebook *stagebook=look->partbooks[partword[j][l][k]][s];
488             if(stagebook){
489               if(decodepart(stagebook,work[j]+i,&vb->opb,
490                             samples_per_partition,0)==-1)goto eopbreak;
491             }
492           }
493         }
494     } 
495   }
496
497  eopbreak:
498   
499   for(j=0;j<ch;j++){
500     for(i=0;i<n;i++)
501       in[j][i]*=work[j][i];
502     for(;i<transend;i++)
503       in[j][i]=0;
504   }
505
506   return(0);
507
508  errout:
509   for(j=0;j<ch;j++)
510     memset(in[j],0,sizeof(float)*transend);
511   return(0);
512 }
513
514 /* residue 0 and 1 are just slight variants of one another. 0 is
515    interleaved, 1 is not */
516 int res0_forward(vorbis_block *vb,vorbis_look_residue *vl,
517             float **in,int ch){
518   return(_01forward(vb,vl,in,ch,_interleaved_testhack,_interleaved_encodepart));
519 }
520
521 int res0_inverse(vorbis_block *vb,vorbis_look_residue *vl,float **in,int ch){
522   return(_01inverse(vb,vl,in,ch,vorbis_book_decodevs));
523 }
524
525 int res1_forward(vorbis_block *vb,vorbis_look_residue *vl,
526                  float **in,int ch){
527   return(_01forward(vb,vl,in,ch,_testhack,_encodepart));
528 }
529
530 int res1_inverse(vorbis_block *vb,vorbis_look_residue *vl,float **in,int ch){
531   return(_01inverse(vb,vl,in,ch,vorbis_book_decodev));
532 }
533
534 vorbis_func_residue residue0_exportbundle={
535   &res0_pack,
536   &res0_unpack,
537   &res0_look,
538   &res0_copy_info,
539   &res0_free_info,
540   &res0_free_look,
541   &res0_forward,
542   &res0_inverse
543 };
544
545 vorbis_func_residue residue1_exportbundle={
546   &res0_pack,
547   &res0_unpack,
548   &res0_look,
549   &res0_copy_info,
550   &res0_free_info,
551   &res0_free_look,
552   &res1_forward,
553   &res1_inverse
554 };