1 /********************************************************************
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. *
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/ *
12 ********************************************************************
14 function: floor backend 0 implementation
15 last mod: $Id: floor0.c,v 1.21 2000/08/19 11:46:28 xiphmont Exp $
17 ********************************************************************/
22 #include "vorbis/codec.h"
27 #include "bookinternal.h"
28 #include "sharedbook.h"
42 vorbis_info_floor0 *vi;
48 /* infrastructure for finding fit */
49 static long _f0_fit(codebook *book,
54 double norm,base=0.,err=0.;
56 double *lsp=workfit+cursor;
58 /* gen a curve for fitting */
59 if(cursor)base=workfit[cursor-1];
60 norm=orig[cursor+dim-1]-base;
64 lsp[i]=(orig[i+cursor]-base);
65 best=_best(book,lsp,1);
67 memcpy(lsp,book->valuelist+best*dim,dim*sizeof(double));
73 /***********************************************/
75 static void free_info(vorbis_info_floor *i){
77 memset(i,0,sizeof(vorbis_info_floor0));
82 static void free_look(vorbis_look_floor *i){
83 vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
85 if(look->linearmap)free(look->linearmap);
86 if(look->lsp_look)free(look->lsp_look);
87 lpc_clear(&look->lpclook);
88 memset(look,0,sizeof(vorbis_look_floor0));
93 static void pack (vorbis_info_floor *i,oggpack_buffer *opb){
94 vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
96 _oggpack_write(opb,info->order,8);
97 _oggpack_write(opb,info->rate,16);
98 _oggpack_write(opb,info->barkmap,16);
99 _oggpack_write(opb,info->ampbits,6);
100 _oggpack_write(opb,info->ampdB,8);
101 _oggpack_write(opb,info->numbooks-1,4);
102 for(j=0;j<info->numbooks;j++)
103 _oggpack_write(opb,info->books[j],8);
106 static vorbis_info_floor *unpack (vorbis_info *vi,oggpack_buffer *opb){
108 vorbis_info_floor0 *info=malloc(sizeof(vorbis_info_floor0));
109 info->order=_oggpack_read(opb,8);
110 info->rate=_oggpack_read(opb,16);
111 info->barkmap=_oggpack_read(opb,16);
112 info->ampbits=_oggpack_read(opb,6);
113 info->ampdB=_oggpack_read(opb,8);
114 info->numbooks=_oggpack_read(opb,4)+1;
116 if(info->order<1)goto err_out;
117 if(info->rate<1)goto err_out;
118 if(info->barkmap<1)goto err_out;
119 if(info->numbooks<1)goto err_out;
121 for(j=0;j<info->numbooks;j++){
122 info->books[j]=_oggpack_read(opb,8);
123 if(info->books[j]<0 || info->books[j]>=vi->books)goto err_out;
131 /* initialize Bark scale and normalization lookups. We could do this
132 with static tables, but Vorbis allows a number of possible
133 combinations, so it's best to do it computationally.
135 The below is authoritative in terms of defining scale mapping.
136 Note that the scale depends on the sampling rate as well as the
137 linear block and mapping sizes */
139 static vorbis_look_floor *look (vorbis_dsp_state *vd,vorbis_info_mode *mi,
140 vorbis_info_floor *i){
143 vorbis_info *vi=vd->vi;
144 vorbis_info_floor0 *info=(vorbis_info_floor0 *)i;
145 vorbis_look_floor0 *look=malloc(sizeof(vorbis_look_floor0));
147 look->n=vi->blocksizes[mi->blockflag]/2;
148 look->ln=info->barkmap;
152 lpc_init(&look->lpclook,look->ln,look->m);
154 /* we choose a scaling constant so that:
155 floor(bark(rate/2-1)*C)=mapped-1
156 floor(bark(rate/2)*C)=mapped */
157 scale=look->ln/toBARK(info->rate/2.);
159 /* the mapping from a linear scale to a smaller bark scale is
160 straightforward. We do *not* make sure that the linear mapping
161 does not skip bark-scale bins; the decoder simply skips them and
162 the encoder may do what it wishes in filling them. They're
163 necessary in some mapping combinations to keep the scale spacing
165 look->linearmap=malloc(look->n*sizeof(int));
166 for(j=0;j<look->n;j++){
167 int val=floor( toBARK((info->rate/2.)/look->n*j)
168 *scale); /* bark numbers represent band edges */
169 if(val>look->ln)val=look->ln; /* guard against the approximation */
170 look->linearmap[j]=val;
173 look->lsp_look=malloc(look->ln*sizeof(double));
174 for(j=0;j<look->ln;j++)
175 look->lsp_look[j]=2*cos(M_PI/look->ln*j);
180 /* less efficient than the decode side (written for clarity). We're
181 not bottlenecked here anyway */
183 double _curve_to_lpc(double *curve,double *lpc,
184 vorbis_look_floor0 *l,long frameno){
185 /* map the input curve to a bark-scale curve for encoding */
188 double *work=alloca(sizeof(double)*mapped);
192 memset(work,0,sizeof(double)*mapped);
194 /* Only the decode side is behavior-specced; for now in the encoder,
195 we select the maximum value of each band as representative (this
196 helps make sure peaks don't go out of range. In error terms,
197 selecting min would make more sense, but the codebook is trained
198 numerically, so we don't actually lose. We'd still want to
199 use the original curve for error and noise estimation */
202 bark=l->linearmap[i];
203 if(work[bark]<curve[i])work[bark]=curve[i];
205 /* If the bark scale is climbing rapidly, some bins may end up
206 going unused. This isn't a waste actually; it keeps the
207 scale resolution even so that the LPC generator has an easy
208 time. However, if we leave the bins empty we lose energy.
209 So, fill 'em in. The decoder does not do anything with he
210 unused bins, so we can fill them anyway we like to end up
211 with a better spectral curve */
213 /* we'll always have a bin zero, so we don't need to guard init */
216 double del=(double)j/span;
217 work[j+last]=work[bark]*del+work[last]*(1.-del);
223 /* If we're over-ranged to avoid edge effects, fill in the end of spectrum gap */
224 for(i=bark+1;i<mapped;i++)
228 { /******************/
233 sprintf(buffer,"Fmask_%d.m",frameno);
234 of=fopen(buffer,"w");
235 for(i=0;i<mapped;i++)
236 fprintf(of,"%g\n",work[i]);
241 return vorbis_lpc_from_curve(work,lpc,&(l->lpclook));
244 /* generate the whole freq response curve of an LPC IIR filter */
246 void _lsp_to_curve(double *curve,double *lsp,double amp,
247 vorbis_look_floor0 *l,char *name,long frameno){
248 /* l->m+1 must be less than l->ln, but guard in case we get a bad stream */
249 double *lcurve=alloca(sizeof(double)*l->ln);
253 memset(curve,0,sizeof(double)*l->n);
256 vorbis_lsp_to_curve(lcurve,l->ln,lsp,l->m,amp,l->lsp_look);
258 for(i=0;i<l->n;i++)curve[i]=lcurve[l->linearmap[i]];
263 static int forward(vorbis_block *vb,vorbis_look_floor *i,
264 double *in,double *out){
266 vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
267 vorbis_info_floor0 *info=look->vi;
268 double *work=alloca((look->ln+look->n)*sizeof(double));
278 sprintf(buffer,"lsp0coeff_%d.vqd",vb->mode);
279 of=fopen(buffer,"a");
282 sprintf(buffer,"lsp0ent_%d.vqd",vb->mode);
283 ef=fopen(buffer,"a");
286 /* our floor comes in on a linear scale; go to a [-Inf...0] dB
287 scale. The curve has to be positive, so we offset it. */
289 for(j=0;j<look->n;j++){
290 double val=todB(in[j])+info->ampdB;
297 /* use 'out' as temp storage */
298 /* Convert our floor to a set of lpc coefficients */
299 amp=sqrt(_curve_to_lpc(work,out,look,seq));
301 /* amp is in the range (0. to ampdB]. Encode that range using
305 long maxval=(1L<<info->ampbits)-1;
307 long val=rint(amp/info->ampdB*maxval);
309 if(val<0)val=0; /* likely */
310 if(val>maxval)val=maxval; /* not bloody likely */
312 _oggpack_write(&vb->opb,val,info->ampbits);
314 amp=(float)val/maxval*info->ampdB;
321 /* the spec supports using one of a number of codebooks. Right
322 now, encode using this lib supports only one */
323 codebook *b=vb->vd->fullbooks+info->books[0];
324 _oggpack_write(&vb->opb,0,_ilog(info->numbooks));
326 /* LSP <-> LPC is orthogonal and LSP quantizes more stably */
327 vorbis_lpc_to_lsp(out,out,look->m);
330 if(vb->W==0){fprintf(stderr,"%d ",seq);}
331 vorbis_lsp_to_lpc(out,work,look->m);
332 _lsp_to_curve(work,work,amp,look,"Ffloor",seq);
333 for(j=0;j<look->n;j++)work[j]-=info->ampdB;
334 _analysis_output("rawfloor",seq,work,look->n,0,0);
337 for(j=0;j<look->m;j++){
342 _analysis_output("rawlsp",seq,work,look->m,0,0);
350 for(j=0;j<look->m;j++){
351 fprintf(of,"%.12g, ",out[j]-last);
360 /* code the spectral envelope, and keep track of the actual
361 quantized values; we don't want creeping error as each block is
362 nailed to the last quantized value of the previous block. */
364 for(j=0;j<look->m;j+=b->dim){
365 int entry=_f0_fit(b,out,work,j);
366 bits+=vorbis_book_encode(b,entry,&vb->opb);
369 fprintf(ef,"%d,\n",entry);
377 for(j=0;j<look->m;j++){
382 _analysis_output("lsp",seq,out,look->m,0,0);
390 /* take the coefficients back to a spectral envelope curve */
391 _lsp_to_curve(out,work,amp,look,"Ffloor",seq++);
392 for(j=0;j<look->n;j++)out[j]= fromdB(out[j]-info->ampdB);
396 memset(out,0,sizeof(double)*look->n);
401 static int inverse(vorbis_block *vb,vorbis_look_floor *i,double *out){
402 vorbis_look_floor0 *look=(vorbis_look_floor0 *)i;
403 vorbis_info_floor0 *info=look->vi;
406 int ampraw=_oggpack_read(&vb->opb,info->ampbits);
407 if(ampraw>0){ /* also handles the -1 out of data case */
408 long maxval=(1<<info->ampbits)-1;
409 double amp=(float)ampraw/maxval*info->ampdB;
410 int booknum=_oggpack_read(&vb->opb,_ilog(info->numbooks));
413 codebook *b=vb->vd->fullbooks+info->books[booknum];
416 memset(out,0,sizeof(double)*look->m);
418 for(j=0;j<look->m;j+=b->dim)
419 if(vorbis_book_decodevs(b,out+j,&vb->opb,1,-1)==-1)goto eop;
421 for(k=0;k<b->dim;k++,j++)out[j]+=last;
425 /* take the coefficients back to a spectral envelope curve */
426 _lsp_to_curve(out,out,amp,look,"",0);
428 for(j=0;j<look->n;j++)out[j]=fromdB(out[j]-info->ampdB);
434 memset(out,0,sizeof(double)*look->n);
439 vorbis_func_floor floor0_exportbundle={
440 &pack,&unpack,&look,&free_info,&free_look,&forward,&inverse