001/* JOrbis
002 * Copyright (C) 2000 ymnk, JCraft,Inc.
003 *  
004 * Written by: 2000 ymnk<ymnk@jcaft.com>
005 *   
006 * Many thanks to 
007 *   Monty <monty@xiph.org> and 
008 *   The XIPHOPHORUS Company http://www.xiph.org/ .
009 * JOrbis has been based on their awesome works, Vorbis codec.
010 *   
011 * This program is free software; you can redistribute it and/or
012 * modify it under the terms of the GNU Library General Public License
013 * as published by the Free Software Foundation; either version 2 of
014 * the License, or (at your option) any later version.
015   
016 * This program is distributed in the hope that it will be useful,
017 * but WITHOUT ANY WARRANTY; without even the implied warranty of
018 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
019 * GNU Library General Public License for more details.
020 * 
021 * You should have received a copy of the GNU Library General Public
022 * License along with this program; if not, write to the Free Software
023 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
024 */
025
026package com.jcraft.jorbis;
027
028import com.jcraft.jogg.*;
029
030class CodeBook{
031  int dim;            // codebook dimensions (elements per vector)
032  int entries;        // codebook entries
033  StaticCodeBook c=new StaticCodeBook();
034
035  float[] valuelist; // list of dim*entries actual entry values
036  int[] codelist;     // list of bitstream codewords for each entry
037  DecodeAux decode_tree;
038
039  // returns the number of bits
040  int encode(int a, Buffer b){
041    b.write(codelist[a], c.lengthlist[a]);
042    return(c.lengthlist[a]);
043  }
044
045  // One the encode side, our vector writers are each designed for a
046  // specific purpose, and the encoder is not flexible without modification:
047  // 
048  // The LSP vector coder uses a single stage nearest-match with no
049  // interleave, so no step and no error return.  This is specced by floor0
050  // and doesn't change.
051  // 
052  // Residue0 encoding interleaves, uses multiple stages, and each stage
053  // peels of a specific amount of resolution from a lattice (thus we want
054  // to match by threshhold, not nearest match).  Residue doesn't *have* to
055  // be encoded that way, but to change it, one will need to add more
056  // infrastructure on the encode side (decode side is specced and simpler)
057
058  // floor0 LSP (single stage, non interleaved, nearest match)
059  // returns entry number and *modifies a* to the quantization value
060  int errorv(float[] a){
061    int best=best(a,1);
062    for(int k=0;k<dim;k++){
063      a[k]=valuelist[best*dim+k];
064    }
065    return(best);
066  }
067
068  // returns the number of bits and *modifies a* to the quantization value
069  int encodev(int best, float[] a, Buffer b){
070    for(int k=0;k<dim;k++){
071      a[k]=valuelist[best*dim+k];
072    }
073    return(encode(best,b));
074  }
075
076  // res0 (multistage, interleave, lattice)
077  // returns the number of bits and *modifies a* to the remainder value
078  int encodevs(float[] a, Buffer b, int step,int addmul){
079    int best=besterror(a,step,addmul);
080    return(encode(best,b));
081  }
082
083  private int[] t=new int[15];  // decodevs_add is synchronized for re-using t.
084  synchronized int decodevs_add(float[]a, int offset, Buffer b, int n){
085    int step=n/dim;
086    int entry;
087    int i,j,o;
088
089    if(t.length<step){
090      t=new int[step];
091    }
092
093    for(i = 0; i < step; i++){
094      entry=decode(b);
095      if(entry==-1)return(-1);
096      t[i]=entry*dim;
097    }
098    for(i=0,o=0;i<dim;i++,o+=step){
099      for(j=0;j<step;j++){
100        a[offset+o+j]+=valuelist[t[j]+i];
101      }
102    }
103
104    return(0);
105  }
106
107  int decodev_add(float[]a, int offset, Buffer b,int n){
108    int i,j,entry;
109    int t;
110
111    if(dim>8){
112      for(i=0;i<n;){
113        entry = decode(b);
114        if(entry==-1)return(-1);
115        t=entry*dim;
116        for(j=0;j<dim;){
117          a[offset+(i++)]+=valuelist[t+(j++)];
118        }
119      }
120    }
121    else{
122      for(i=0;i<n;){
123        entry=decode(b);
124        if(entry==-1)return(-1);
125        t=entry*dim;
126        j=0;
127        switch(dim){
128        case 8:
129          a[offset+(i++)]+=valuelist[t+(j++)];
130        case 7:
131          a[offset+(i++)]+=valuelist[t+(j++)];
132        case 6:
133          a[offset+(i++)]+=valuelist[t+(j++)];
134        case 5:
135          a[offset+(i++)]+=valuelist[t+(j++)];
136        case 4:
137          a[offset+(i++)]+=valuelist[t+(j++)];
138        case 3:
139          a[offset+(i++)]+=valuelist[t+(j++)];
140        case 2:
141          a[offset+(i++)]+=valuelist[t+(j++)];
142        case 1:
143          a[offset+(i++)]+=valuelist[t+(j++)];
144        case 0:
145        break;
146        }
147      }
148    }    
149    return(0);
150  }
151
152  int decodev_set(float[] a,int offset, Buffer b, int n){
153    int i,j,entry;
154    int t;
155
156    for(i=0;i<n;){
157      entry = decode(b);
158      if(entry==-1)return(-1);
159      t=entry*dim;
160      for(j=0;j<dim;){
161        a[offset+i++]=valuelist[t+(j++)];
162      }
163    }
164    return(0);
165  }
166
167  int decodevv_add(float[][] a, int offset,int ch, Buffer b,int n){
168    int i,j,k,entry;
169    int chptr=0;
170    //System.out.println("decodevv_add: a="+a+",b="+b+",valuelist="+valuelist);
171
172    for(i=offset/ch;i<(offset+n)/ch;){
173      entry = decode(b);
174      if(entry==-1)return(-1);
175
176      int t = entry*dim;
177      for(j=0;j<dim;j++){
178        a[chptr++][i]+=valuelist[t+j];
179        if(chptr==ch){
180          chptr=0;
181          i++;
182        }
183      }
184    }
185    return(0);
186  }
187
188
189  // Decode side is specced and easier, because we don't need to find
190  // matches using different criteria; we simply read and map.  There are
191  // two things we need to do 'depending':
192  //   
193  // We may need to support interleave.  We don't really, but it's
194  // convenient to do it here rather than rebuild the vector later.
195  //
196  // Cascades may be additive or multiplicitive; this is not inherent in
197  // the codebook, but set in the code using the codebook.  Like
198  // interleaving, it's easiest to do it here.  
199  // stage==0 -> declarative (set the value)
200  // stage==1 -> additive
201  // stage==2 -> multiplicitive
202
203  // returns the entry number or -1 on eof
204  int decode(Buffer b){
205    int ptr=0;
206    DecodeAux t=decode_tree;
207    int lok=b.look(t.tabn);
208    //System.err.println(this+" "+t+" lok="+lok+", tabn="+t.tabn);
209
210    if(lok>=0){
211      ptr=t.tab[lok];
212      b.adv(t.tabl[lok]);
213      if(ptr<=0){
214        return -ptr;
215      }
216    }
217    do{
218      switch(b.read1()){
219      case 0:
220        ptr=t.ptr0[ptr];
221        break;
222      case 1:
223        ptr=t.ptr1[ptr];
224        break;
225      case -1:
226      default:
227        return(-1);
228      }
229    }
230    while(ptr>0);
231    return(-ptr);
232  }
233
234  // returns the entry number or -1 on eof
235  int decodevs(float[] a, int index, Buffer b, int step,int addmul){
236    int entry=decode(b);
237    if(entry==-1)return(-1);
238    switch(addmul){
239    case -1:
240      for(int i=0,o=0;i<dim;i++,o+=step)
241        a[index+o]=valuelist[entry*dim+i];
242      break;
243    case 0:
244      for(int i=0,o=0;i<dim;i++,o+=step)
245        a[index+o]+=valuelist[entry*dim+i];
246      break;
247    case 1:
248      for(int i=0,o=0;i<dim;i++,o+=step)
249        a[index+o]*=valuelist[entry*dim+i];
250      break;
251    default:
252      System.err.println("CodeBook.decodeves: addmul="+addmul); 
253    }
254    return(entry);
255  }
256
257  int best(float[] a, int step){
258    EncodeAuxNearestMatch nt=c.nearest_tree;
259    EncodeAuxThreshMatch tt=c.thresh_tree;
260    int ptr=0;
261
262    // we assume for now that a thresh tree is the only other possibility
263    if(tt!=null){
264      int index=0;
265      // find the quant val of each scalar
266      for(int k=0,o=step*(dim-1);k<dim;k++,o-=step){
267        int i;
268        // linear search the quant list for now; it's small and although
269        // with > 8 entries, it would be faster to bisect, this would be
270        // a misplaced optimization for now
271        for(i=0;i<tt.threshvals-1;i++){
272          if(a[o]<tt.quantthresh[i]){
273            break;
274          }
275        }
276        index=(index*tt.quantvals)+tt.quantmap[i];
277      }
278      // regular lattices are easy :-)
279      if(c.lengthlist[index]>0){
280        // is this unused?  If so, we'll
281        // use a decision tree after all
282        // and fall through
283        return(index);
284      }
285    }
286    if(nt!=null){
287      // optimized using the decision tree
288      while(true){
289        float c=0.f;
290        int p=nt.p[ptr];
291        int q=nt.q[ptr];
292        for(int k=0,o=0;k<dim;k++,o+=step){
293          c+=(valuelist[p+k]-valuelist[q+k])*
294             (a[o]-(valuelist[p+k]+valuelist[q+k])*.5);
295        }
296        if(c>0.){ // in A
297          ptr= -nt.ptr0[ptr];
298        }
299        else{     // in B
300          ptr= -nt.ptr1[ptr];
301        }
302        if(ptr<=0)break;
303      }
304      return(-ptr);
305    }
306
307    // brute force it!
308    {
309      int besti=-1;
310      float best=0.f;
311      int e=0;
312      for(int i=0;i<entries;i++){
313        if(c.lengthlist[i]>0){
314          float _this=dist(dim, valuelist, e, a, step);
315          if(besti==-1 || _this<best){
316            best=_this;
317            besti=i;
318          }
319        }
320        e+=dim;
321      }
322      return(besti);
323    }
324  }
325
326  // returns the entry number and *modifies a* to the remainder value
327  int besterror(float[] a, int step, int addmul){
328    int best=best(a,step);
329    switch(addmul){
330    case 0:
331      for(int i=0,o=0;i<dim;i++,o+=step)
332        a[o]-=valuelist[best*dim+i];
333      break;
334    case 1:
335      for(int i=0,o=0;i<dim;i++,o+=step){
336        float val=valuelist[best*dim+i];
337        if(val==0){
338          a[o]=0;
339        }else{
340          a[o]/=val;
341        }
342      }
343      break;
344    }
345    return(best);
346  }
347
348  void clear(){
349    // static book is not cleared; we're likely called on the lookup and
350    // the static codebook belongs to the info struct
351    //if(decode_tree!=null){
352    //  free(b->decode_tree->ptr0);
353    //  free(b->decode_tree->ptr1);
354    //  memset(b->decode_tree,0,sizeof(decode_aux));
355    //  free(b->decode_tree);
356    //}
357    //if(valuelist!=null)free(b->valuelist);
358    //if(codelist!=null)free(b->codelist);
359    //memset(b,0,sizeof(codebook));
360  }
361
362  private static float dist(int el, float[] ref, int index, float[] b, int step){
363    float acc=(float)0.;
364    for(int i=0; i<el; i++){
365      float val=(ref[index+i]-b[i*step]);
366      acc+=val*val;
367    }
368    return(acc);
369  }
370
371/*
372  int init_encode(StaticCodeBook s){
373    //memset(c,0,sizeof(codebook));
374    c=s;
375    entries=s.entries;
376    dim=s.dim;
377    codelist=make_words(s.lengthlist, s.entries);
378    valuelist=s.unquantize();
379    return(0);
380  }
381*/
382
383  int init_decode(StaticCodeBook s){
384    //memset(c,0,sizeof(codebook));
385    c=s;
386    entries=s.entries;
387    dim=s.dim;
388    valuelist=s.unquantize();
389
390    decode_tree=make_decode_tree();
391    if(decode_tree==null){
392      //goto err_out;
393      clear();
394      return(-1);
395    }
396    return(0);
397//  err_out:
398//    vorbis_book_clear(c);
399//    return(-1);
400  }
401
402  // given a list of word lengths, generate a list of codewords.  Works
403  // for length ordered or unordered, always assigns the lowest valued
404  // codewords first.  Extended to handle unused entries (length 0)
405  static int[] make_words(int[] l, int n){
406    int[] marker=new int[33];
407    int[] r=new int[n];
408    //memset(marker,0,sizeof(marker));
409
410    for(int i=0;i<n;i++){
411      int length=l[i];
412      if(length>0){
413        int entry=marker[length];
414      
415        // when we claim a node for an entry, we also claim the nodes
416        // below it (pruning off the imagined tree that may have dangled
417        // from it) as well as blocking the use of any nodes directly
418        // above for leaves
419      
420        // update ourself
421        if(length<32 && (entry>>>length)!=0){
422          // error condition; the lengths must specify an overpopulated tree
423          //free(r);
424          return(null);
425        }
426        r[i]=entry;
427    
428        // Look to see if the next shorter marker points to the node
429        // above. if so, update it and repeat.
430        {
431          for(int j=length;j>0;j--){
432            if((marker[j]&1)!=0){
433              // have to jump branches
434              if(j==1)marker[1]++;
435              else marker[j]=marker[j-1]<<1;
436              break; // invariant says next upper marker would already
437                     // have been moved if it was on the same path
438            }
439            marker[j]++;
440          }
441        }
442      
443        // prune the tree; the implicit invariant says all the longer
444        // markers were dangling from our just-taken node.  Dangle them
445        // from our *new* node.
446        for(int j=length+1;j<33;j++){
447          if((marker[j]>>>1) == entry){
448            entry=marker[j];
449            marker[j]=marker[j-1]<<1;
450          }
451          else{
452            break;
453          }
454        }    
455      }
456    }
457
458    // bitreverse the words because our bitwise packer/unpacker is LSb
459    // endian
460    for(int i=0;i<n;i++){
461      int temp=0;
462      for(int j=0;j<l[i];j++){
463        temp<<=1;
464        temp|=(r[i]>>>j)&1;
465      }
466      r[i]=temp;
467    }
468
469    return(r);
470  }
471
472  // build the decode helper tree from the codewords 
473  DecodeAux make_decode_tree(){
474    int top=0;
475    DecodeAux t=new DecodeAux();
476    int[] ptr0=t.ptr0=new int[entries*2];
477    int[] ptr1=t.ptr1=new int[entries*2];
478    int[] codelist=make_words(c.lengthlist, c.entries);
479
480    if(codelist==null)return(null);
481    t.aux=entries*2;
482
483    for(int i=0;i<entries;i++){
484      if(c.lengthlist[i]>0){
485        int ptr=0;
486        int j;
487        for(j=0;j<c.lengthlist[i]-1;j++){
488          int bit=(codelist[i]>>>j)&1;
489          if(bit==0){
490            if(ptr0[ptr]==0){
491              ptr0[ptr]=++top;
492            }
493            ptr=ptr0[ptr];
494          }
495          else{
496            if(ptr1[ptr]==0){
497              ptr1[ptr]= ++top;
498            }
499            ptr=ptr1[ptr];
500          }
501        }
502
503        if(((codelist[i]>>>j)&1)==0){ ptr0[ptr]=-i; }
504        else{ ptr1[ptr]=-i; }
505
506      }
507    }
508    //free(codelist);
509
510    t.tabn = ilog(entries)-4;
511
512    if(t.tabn<5)t.tabn=5;
513    int n = 1<<t.tabn;
514    t.tab = new int[n];
515    t.tabl = new int[n];
516    for(int i = 0; i < n; i++){
517      int p = 0;
518      int j=0;
519      for(j = 0; j < t.tabn && (p > 0 || j == 0); j++){
520        if ((i&(1<<j))!=0){
521          p = ptr1[p];
522        }
523        else{
524          p = ptr0[p];
525        }
526      }
527      t.tab[i]=p;  // -code
528      t.tabl[i]=j; // length 
529    }
530
531    return(t);
532  }
533
534  private static int ilog(int v){
535    int ret=0;
536    while(v!=0){
537      ret++;
538      v>>>=1;
539    }
540    return(ret);
541  }
542
543/*
544  // TEST
545  // Simple enough; pack a few candidate codebooks, unpack them.  Code a
546  // number of vectors through (keeping track of the quantized values),
547  // and decode using the unpacked book.  quantized version of in should
548  // exactly equal out
549
550  //#include "vorbis/book/lsp20_0.vqh"
551  //#include "vorbis/book/lsp32_0.vqh"
552  //#include "vorbis/book/res0_1a.vqh"
553  static final int TESTSIZE=40;
554
555  static float[] test1={
556    0.105939,
557    0.215373,
558    0.429117,
559    0.587974,
560
561    0.181173,
562    0.296583,
563    0.515707,
564    0.715261,
565
566    0.162327,
567    0.263834,
568    0.342876,
569    0.406025,
570
571    0.103571,
572    0.223561,
573    0.368513,
574    0.540313,
575
576    0.136672,
577    0.395882,
578    0.587183,
579    0.652476,
580
581    0.114338,
582    0.417300,
583    0.525486,
584    0.698679,
585
586    0.147492,
587    0.324481,
588    0.643089,
589    0.757582,
590
591    0.139556,
592    0.215795,
593    0.324559,
594    0.399387,
595
596    0.120236,
597    0.267420,
598    0.446940,
599    0.608760,
600
601    0.115587,
602    0.287234,
603    0.571081,
604    0.708603,
605  };
606
607  static float[] test2={
608    0.088654,
609    0.165742,
610    0.279013,
611    0.395894,
612
613    0.110812,
614    0.218422,
615    0.283423,
616    0.371719,
617
618    0.136985,
619    0.186066,
620    0.309814,
621    0.381521,
622
623    0.123925,
624    0.211707,
625    0.314771,
626    0.433026,
627
628    0.088619,
629    0.192276,
630    0.277568,
631    0.343509,
632
633    0.068400,
634    0.132901,
635    0.223999,
636    0.302538,
637
638    0.202159,
639    0.306131,
640    0.360362,
641    0.416066,
642
643    0.072591,
644    0.178019,
645    0.304315,
646    0.376516,
647
648    0.094336,
649    0.188401,
650    0.325119,
651    0.390264,
652
653    0.091636,
654    0.223099,
655    0.282899,
656    0.375124,
657  };
658
659  static float[] test3={
660    0,1,-2,3,4,-5,6,7,8,9,
661    8,-2,7,-1,4,6,8,3,1,-9,
662    10,11,12,13,14,15,26,17,18,19,
663    30,-25,-30,-1,-5,-32,4,3,-2,0};
664
665//  static_codebook *testlist[]={&_vq_book_lsp20_0,
666//                             &_vq_book_lsp32_0,
667//                             &_vq_book_res0_1a,NULL};
668  static[][] float testvec={test1,test2,test3};
669
670  static void main(String[] arg){
671    Buffer write=new Buffer();
672    Buffer read=new Buffer();
673    int ptr=0;
674    write.writeinit();
675  
676    System.err.println("Testing codebook abstraction...:");
677
678    while(testlist[ptr]!=null){
679      CodeBook c=new CodeBook();
680      StaticCodeBook s=new StaticCodeBook();;
681      float *qv=alloca(sizeof(float)*TESTSIZE);
682      float *iv=alloca(sizeof(float)*TESTSIZE);
683      memcpy(qv,testvec[ptr],sizeof(float)*TESTSIZE);
684      memset(iv,0,sizeof(float)*TESTSIZE);
685
686      System.err.print("\tpacking/coding "+ptr+"... ");
687
688      // pack the codebook, write the testvector
689      write.reset();
690      vorbis_book_init_encode(&c,testlist[ptr]); // get it into memory
691                                                 // we can write
692      vorbis_staticbook_pack(testlist[ptr],&write);
693      System.err.print("Codebook size "+write.bytes()+" bytes... ");
694      for(int i=0;i<TESTSIZE;i+=c.dim){
695        vorbis_book_encodev(&c,qv+i,&write);
696      }
697      c.clear();
698    
699      System.err.print("OK.\n");
700      System.err.print("\tunpacking/decoding "+ptr+"... ");
701
702      // transfer the write data to a read buffer and unpack/read
703      _oggpack_readinit(&read,_oggpack_buffer(&write),_oggpack_bytes(&write));
704      if(s.unpack(read)){
705        System.err.print("Error unpacking codebook.\n");
706        System.exit(1);
707      }
708      if(vorbis_book_init_decode(&c,&s)){
709        System.err.print("Error initializing codebook.\n");
710        System.exit(1);
711      }
712      for(int i=0;i<TESTSIZE;i+=c.dim){
713        if(vorbis_book_decodevs(&c,iv+i,&read,1,-1)==-1){
714          System.err.print("Error reading codebook test data (EOP).\n");
715          System.exit(1);
716        }
717      }
718      for(int i=0;i<TESTSIZE;i++){
719        if(fabs(qv[i]-iv[i])>.000001){
720          System.err.print("read ("+iv[i]+") != written ("+qv[i]+") at position ("+i+")\n");
721          System.exit(1);
722        }
723      }
724
725      System.err.print("OK\n");
726      ptr++;
727    }
728    // The above is the trivial stuff; 
729    // now try unquantizing a log scale codebook
730  }
731*/
732}
733
734class DecodeAux{
735  int[] tab;
736  int[] tabl;
737  int tabn;
738
739  int[] ptr0;
740  int[] ptr1;
741  int   aux;        // number of tree entries
742}