- if(i%2)lpc[j]+=lpc[j]*r;
-
- error*=1.0-r*r;
- }
-
- /* we need the error value to know how big an impulse to hit the
- filter with later */
-
- return error;
-}
-
-/* Input : n element envelope spectral curve
- Output: m lpc coefficients, excitation energy */
-
-double vorbis_lpc_from_spectrum(double *curve,double *lpc,lpc_lookup *l){
- int n=l->ln;
- int m=l->m;
- double *work=alloca(sizeof(double)*(n+n));
- double fscale=.5/n;
- int i,j;
-
- /* input is a real curve. make it complex-real */
- /* This mixes phase, but the LPC generation doesn't care. */
- for(i=0;i<n;i++){
- work[i*2]=curve[i]*fscale;
- work[i*2+1]=0;
- }
-
- n*=2;
- drft_backward(&l->fft,work);
-
- /* The autocorrelation will not be circular. Shift, else we lose
- most of the power in the edges. */
-
- for(i=0,j=n/2;i<n/2;){
- double temp=work[i];
- work[i++]=work[j];
- work[j++]=temp;
- }
-
- return(vorbis_lpc_from_data(work,lpc,n,m));
-}
-
-/* initialize Bark scale and normalization lookups. We could do this
- with static tables, but Vorbis allows a number of possible
- combinations, so it's best to do it computationally.
-
- The below is authoritative in terms of defining scale mapping.
- Note that the scale depends on the sampling rate as well as the
- linear block and mapping sizes (note that for a given sample rate
- and block size, there's generally a fairly obviously optimal
- mapping size */
-
-void lpc_init(lpc_lookup *l,int n, long mapped, long rate, int m){
- int i;
- double scale;
- memset(l,0,sizeof(lpc_lookup));
-
- l->n=n;
- l->ln=mapped;
- l->m=m;
-
- l->linearmap=malloc(n*sizeof(int));
- l->barknorm=malloc(mapped*sizeof(double));
-
- /* we choose a scaling constant so that:
- floor(bark(rate-1)*C)=mapped-1
- floor(bark(rate)*C)=mapped */
-
- scale=mapped/fBARK(rate);
-
- /* the mapping from a linear scale to a smaller bark scale is
- straightforward with a single catch; make sure not to skip any
- bark-scale bins. In order to do this, we assign map_N = min
- (map_N-1 + 1, bark(N)) */
- {
- int last=-1;
- for(i=0;i<n;i++){
- int val=floor( fBARK(((double)rate)/n*i) *scale); /* bark numbers
- represent
- band edges */
- if(val>=mapped)val=mapped; /* guard against the approximation */
- if(val>last+1)val=last+1;
- l->linearmap[i]=val;
- last=val;
- }
- }
-
- /* 'Normalization' is just making sure that power isn't lost in the
- log scale by virtue of compressing the scale in higher
- frequencies. We figure the weight of bands in proportion to
- their linear/bark width ratio below, again, authoritatively. We
- use computed width (not the number of actual bins above) for
- smoothness in the scale; they should agree closely unless the
- encoder chose parameters poorly (and got a bark scale that would
- have had lots of skipped bins) */
-
- for(i=0;i<mapped;i++)
- l->barknorm[i]=iBARK((i+1)/scale)-iBARK(i/scale);
-
- /* we cheat decoding the LPC spectrum via FFTs */
-
- drft_init(&l->fft,mapped*2);
-
-}
-
-void lpc_clear(lpc_lookup *l){
- if(l){
- if(l->barknorm)free(l->barknorm);
- if(l->linearmap)free(l->linearmap);
- drft_clear(&l->fft);
- }
-}
-
-
-/* less efficient than the decode side (written for clarity). We're
- not bottlenecked here anyway */
-
-double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l){
- /* map the input curve to a bark-scale curve for encoding */
-
- int mapped=l->ln;
- double *work=alloca(sizeof(double)*mapped);
- int i;
-
- memset(work,0,sizeof(double)*mapped);
-
- /* Only the decode side is behavior-specced; for now in the encoder,
- we select the maximum value of each band as representative (this
- helps make sure peaks don't go out of range. In error terms,
- selecting min would make more sense, but the codebook is trained
- numerically, so we don't lose in encoding. We'd still want to
- use the original curve for error and noise estimation */
-
- for(i=0;i<l->n;i++){
- int bark=l->linearmap[i];
- if(work[bark]<curve[i])work[bark]=curve[i];
- }
- for(i=0;i<mapped;i++)work[i]*=l->barknorm[i];
-
-#ifdef ANALYSIS
- {
- int j;
- FILE *out;
- char buffer[80];
- static int frameno=0;
-
- sprintf(buffer,"prelpc%d.m",frameno);
- out=fopen(buffer,"w+");
- for(j=0;j<l->n;j++)
- fprintf(out,"%g\n",curve[j]);
- fclose(out);
- sprintf(buffer,"preloglpc%d.m",frameno++);
- out=fopen(buffer,"w+");
- for(j=0;j<l->ln;j++)
- fprintf(out,"%g\n",work[j]);
- fclose(out);
- }
-#endif
-
- return vorbis_lpc_from_spectrum(work,lpc,l);
-}
-
-
-/* One can do this the long way by generating the transfer function in
- the time domain and taking the forward FFT of the result. The
- results from direct calculation are cleaner and faster.