Time domain clamping backed out.
Monty
svn path=/trunk/vorbis/; revision=153
function: single-block PCM analysis
author: Monty <xiphmont@mit.edu>
modifications by: Monty
- last modification date: Oct 15 1999
+ last modification date: Oct 21 1999
********************************************************************/
_oggpack_reset(opb);
/* Encode the packet type */
_oggpack_write(opb,0,1);
+
/* Encode the block size */
_oggpack_write(opb,vb->W,1);
+ if(vb->W){
+ _oggpack_write(opb,vb->lW,1);
+ _oggpack_write(opb,vb->nW,1);
+ }
- /* we have the preecho metrics; decide what to do with them */
- _ve_envelope_sparsify(vb);
- _ve_envelope_apply(vb,0);
-
- /* Encode the envelope */
- if(_ve_envelope_encode(vb))return(-1);
+ /* No envelope encoding yet */
+ _oggpack_write(opb,0,1);
/* time domain PCM -> MDCT domain */
for(i=0;i<vi->channels;i++)
FILE *out;
char buffer[80];
- sprintf(buffer,"spectrum.m");
+ sprintf(buffer,"Aspectrum%d.m",vb->sequence);
out=fopen(buffer,"w+");
for(j=0;j<n/2;j++)
fprintf(out,"%g\n",vb->pcm[i][j]);
fclose(out);
- sprintf(buffer,"noise.m");
+ sprintf(buffer,"Anoise%d.m",vb->sequence);
out=fopen(buffer,"w+");
for(j=0;j<n/2;j++)
fprintf(out,"%g\n",floor[j]);
FILE *out;
char buffer[80];
- sprintf(buffer,"premask.m");
+ sprintf(buffer,"Apremask%d.m",vb->sequence);
out=fopen(buffer,"w+");
for(j=0;j<n/2;j++)
fprintf(out,"%g\n",floor[j]);
FILE *out;
char buffer[80];
- sprintf(buffer,"mask.m");
+ sprintf(buffer,"Alpc%d.m",vb->sequence);
+ out=fopen(buffer,"w+");
+ for(j=0;j<vl->m;j++)
+ fprintf(out,"%g\n",lpc[j]);
+ fclose(out);
+
+ sprintf(buffer,"Alsp%d.m",vb->sequence);
+ out=fopen(buffer,"w+");
+ for(j=0;j<vl->m;j++)
+ fprintf(out,"%g\n",lsp[j]);
+ fclose(out);
+
+ sprintf(buffer,"Amask%d.m",vb->sequence);
out=fopen(buffer,"w+");
for(j=0;j<n/2;j++)
fprintf(out,"%g\n",curve[j]);
fclose(out);
- sprintf(buffer,"res.m");
+ sprintf(buffer,"Ares%d.m",vb->sequence);
out=fopen(buffer,"w+");
for(j=0;j<n/2;j++)
fprintf(out,"%g\n",vb->pcm[i][j]);
op->b_o_s=0;
op->e_o_s=vb->eofflag;
op->frameno=vb->frameno;
+ op->packetno=vb->sequence; /* for sake of completeness */
return(0);
}
function: PCM data vector blocking, windowing and dis/reassembly
author: Monty <xiphmont@mit.edu>
modifications by: Monty
- last modification date: Oct 15 1999
+ last modification date: Oct 22 1999
Handle windowing, overlap-add, etc of the PCM vectors. This is made
more amusing by Vorbis' current two allowed block sizes.
********************************************************************/
+#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "codec.h"
vb->pcm_storage=vi->blocksize[1];
vb->pcm_channels=vi->channels;
- vb->mult_storage=vi->blocksize[1]/vi->envelopesa;
- vb->mult_channels=vi->envelopech;
vb->floor_channels=vi->floorch;
vb->floor_storage=max(vi->floororder[0],vi->floororder[1]);
for(i=0;i<vb->pcm_channels;i++)
vb->pcm[i]=malloc(vb->pcm_storage*sizeof(double));
- vb->mult=malloc(vb->mult_channels*sizeof(double *));
- for(i=0;i<vb->mult_channels;i++)
- vb->mult[i]=malloc(vb->mult_storage*sizeof(double));
-
vb->lsp=malloc(vb->floor_channels*sizeof(double *));
vb->lpc=malloc(vb->floor_channels*sizeof(double *));
vb->amp=malloc(vb->floor_channels*sizeof(double));
free(vb->pcm[i]);
free(vb->pcm);
}
- if(vb->mult){
- for(i=0;i<vb->mult_channels;i++)
- free(vb->mult[i]);
- free(vb->mult);
- }
if(vb->vd->analysisp)
_oggpack_writeclear(&vb->opb);
memset(v,0,sizeof(vorbis_dsp_state));
v->vi=vi;
- _ve_envelope_init(&v->ve,vi->envelopesa);
mdct_init(&v->vm[0],vi->blocksize[0]);
mdct_init(&v->vm[1],vi->blocksize[1]);
v->pcm[i]=calloc(v->pcm_storage,sizeof(double));
}
- /* Initialize the envelope multiplier storage */
-
- if(vi->envelopech){
- v->envelope_storage=v->pcm_storage/vi->envelopesa;
- v->multipliers=calloc(vi->envelopech,sizeof(double *));
- {
- int i;
- for(i=0;i<vi->envelopech;i++){
- v->multipliers[i]=calloc(v->envelope_storage,sizeof(double));
- }
- }
- }
-
/* all 1 (large block) or 0 (small block) */
/* explicitly set for the sake of clarity */
v->lW=0; /* previous window size */
v->centerW=vi->blocksize[1]/2;
v->pcm_current=v->centerW;
- v->envelope_current=v->centerW/vi->envelopesa;
return(0);
}
int vorbis_analysis_init(vorbis_dsp_state *v,vorbis_info *vi){
_vds_shared_init(v,vi);
+ /* Initialize the envelope multiplier storage */
+
+ v->envelope_storage=v->pcm_storage/vi->envelopesa;
+ v->multipliers=calloc(v->envelope_storage,sizeof(double));
+ _ve_envelope_init(&v->ve,vi->envelopesa);
+
/* the coder init is different for read/write */
v->analysisp=1;
_vp_psy_init(&v->vp[0],vi,vi->blocksize[0]/2);
lpc_init(&v->vbal[1],vi->blocksize[1]/2,256,
vi->balanceorder,vi->balanceoctaves,1);*/
+ v->envelope_current=v->centerW/vi->envelopesa;
+
return(0);
}
free(v->pcm);
if(v->pcmret)free(v->pcmret);
}
- if(v->multipliers){
- for(i=0;i<vi->envelopech;i++)
- if(v->multipliers[i])free(v->multipliers[i]);
- free(v->multipliers);
- }
+ if(v->multipliers)free(v->multipliers);
+
_ve_envelope_clear(&v->ve);
mdct_clear(&v->vm[0]);
mdct_clear(&v->vm[1]);
for(i=0;i<vi->channels;i++){
v->pcm[i]=realloc(v->pcm[i],v->pcm_storage*sizeof(double));
}
- for(i=0;i<vi->envelopech;i++){
- v->multipliers[i]=realloc(v->multipliers[i],
- v->envelope_storage*sizeof(double));
- }
+ v->multipliers=realloc(v->multipliers,v->envelope_storage*sizeof(double));
}
for(i=0;i<vi->channels;i++)
/* do the deltas, envelope shaping, pre-echo and determine the size of
the next block on which to continue analysis */
int vorbis_analysis_blockout(vorbis_dsp_state *v,vorbis_block *vb){
- int i,j;
+ int i;
vorbis_info *vi=v->vi;
long beginW=v->centerW-vi->blocksize[v->W]/2,centerNext;
- long beginM=beginW/vi->envelopesa;
/* check to see if we're done... */
if(v->eofflag==-1)return(0);
if(v->pcm_current/vi->envelopesa>v->envelope_current){
/* This generates the multipliers, but does not sparsify the vector.
That's done by block before coding */
- _ve_envelope_multipliers(v);
+ _ve_envelope_deltas(v);
}
/* By our invariant, we have lW, W and centerW set. Search for
the next boundary so we can determine nW (the next window size)
which lets us compute the shape of the current block's window */
-
- /* overconserve for now; any block with a non-placeholder multiplier
- should be minimal size. We can be greedy and only look at nW size */
if(vi->blocksize[0]<vi->blocksize[1]){
i=v->centerW;
i/=vi->envelopesa;
- for(;i<v->envelope_current;i++){
- for(j=0;j<vi->envelopech;j++)
- if(v->multipliers[j][i-1]*vi->preecho_thresh<
- v->multipliers[j][i])break;
- if(j<vi->envelopech)break;
+ for(;i<v->envelope_current-1;i++){
+ /* Compare last with current; do we have an abrupt energy change? */
+
+ if(v->multipliers[i-1]*vi->preecho_thresh<
+ v->multipliers[i])break;
+
+ /* because the overlapping nature of the delta finding
+ 'smears' the energy cliffs, also compare completely
+ unoverlapped areas just in case the plosive happened in an
+ unlucky place */
+
+ if(v->multipliers[i-1]*vi->preecho_thresh<
+ v->multipliers[i+1])break;
+
}
- if(i<v->envelope_current){
+ if(i<v->envelope_current-1){
/* Ooo, we hit a multiplier. Is it beyond the boundary to make the
upcoming block large ? */
long largebound;
if(v->W)
- largebound=v->centerW+vi->blocksize[1];
+ /* min boundary; nW large, next small */
+ largebound=v->centerW+vi->blocksize[1]*3/4+vi->blocksize[0]/4;
else
- largebound=v->centerW+vi->blocksize[0]/4+vi->blocksize[1]*3/4;
+ /* min boundary; nW large, next small */
+ largebound=v->centerW+vi->blocksize[0]/2+vi->blocksize[1]/2;
largebound/=vi->envelopesa;
if(i>=largebound)
v->nW=1;
}
}else
- v->nW=1;
+ v->nW=0;
/* Do we actually have enough data *now* for the next block? The
reason to check is that if we had no multipliers, that could
- simply been due to running out of data. In that case, we don;t
+ simply been due to running out of data. In that case, we don't
know the size of the next block for sure and we need that now to
figure out the window shape of this block */
centerNext=v->centerW+vi->blocksize[v->W]/4+vi->blocksize[v->nW]/4;
{
+ /* center of next block + next block maximum right side. Note
+ that the next block needs an additional vi->envelopesa samples
+ to actually be written (for the last multiplier), but we didn't
+ need that to determine its size */
+
long blockbound=centerNext+vi->blocksize[v->nW]/2;
if(v->pcm_current<blockbound)return(0); /* not enough data yet */
}
-
- /* fill in the block */
- vb->lW=v->lW;
- vb->W=v->W;
- vb->nW=v->nW;
+
+ /* fill in the block. Note that for a short window, lW and nW are *short*
+ regardless of actual settings in the stream */
+ fprintf(stderr,"%d",v->W);
+ if(v->W){
+ vb->lW=v->lW;
+ vb->W=v->W;
+ vb->nW=v->nW;
+ }else{
+ vb->lW=0;
+ vb->W=v->W;
+ vb->nW=0;
+ }
vb->vd=v;
-
- vb->pcmend=vi->blocksize[v->W];
- vb->multend=vb->pcmend / vi->envelopesa;
-
+ vb->sequence=v->sequence;
+ vb->frameno=v->frameno;
+
/* copy the vectors */
+ vb->pcmend=vi->blocksize[v->W];
for(i=0;i<vi->channels;i++)
memcpy(vb->pcm[i],v->pcm[i]+beginW,vi->blocksize[v->W]*sizeof(double));
- for(i=0;i<vi->envelopech;i++)
- memcpy(vb->mult[i],v->multipliers[i]+beginM,vi->blocksize[v->W]/
- vi->envelopesa*sizeof(double));
-
- vb->sequence=v->sequence;
- vb->frameno=v->frameno;
-
+
/* handle eof detection: eof==0 means that we've not yet received EOF
eof>0 marks the last 'real' sample in pcm[]
eof<0 'no more to do'; doesn't get here */
must be multiples of samples_per_envelope_step (minimum
multiple is 2) */
+ v->pcm_current-=movementW;
+ v->envelope_current-=movementM;
+
for(i=0;i<vi->channels;i++)
memmove(v->pcm[i],v->pcm[i]+movementW,
- (v->pcm_current-movementW)*sizeof(double));
+ v->pcm_current*sizeof(double));
- for(i=0;i<vi->envelopech;i++){
- memmove(v->multipliers[i],v->multipliers[i]+movementM,
- (v->envelope_current-movementM)*sizeof(double));
- }
-
- v->pcm_current-=movementW;
- v->envelope_current-=movementM;
+ memmove(v->multipliers,v->multipliers+movementM,
+ v->envelope_current*sizeof(double));
v->lW=v->W;
v->W=v->nW;
}
int vorbis_synthesis_init(vorbis_dsp_state *v,vorbis_info *vi){
- int temp=vi->envelopech;
- vi->envelopech=0; /* we don't need multiplier buffering in syn */
_vds_shared_init(v,vi);
- vi->envelopech=temp;
/* Yes, wasteful to have four lookups. This will get collapsed once
things crystallize */
lpc_init(&v->vl[0],vi->blocksize[0]/2,vi->blocksize[0]/2,
vi->floororder[0],vi->flooroctaves[0],0);
- lpc_init(&v->vl[1],vi->blocksize[1]/2,vi->blocksize[0]/2,
+ lpc_init(&v->vl[1],vi->blocksize[1]/2,vi->blocksize[1]/2,
vi->floororder[1],vi->flooroctaves[1],0);
/*lpc_init(&v->vbal[0],vi->blocksize[0]/2,256,
vi->balanceorder,vi->balanceoctaves,0);
return(0);
}
-/* Unike in analysis, the window is only partially applied. Envelope
- is previously applied (the whole envelope, if any, is shipped in
- each block) */
+/* Unike in analysis, the window is only partially applied for each
+ block. The time domain envelope is not yet handled at the point of
+ calling (as it relies on the previous block). */
int vorbis_synthesis_blockin(vorbis_dsp_state *v,vorbis_block *vb){
vorbis_info *vi=v->vi;
- /* Shift out any PCM that we returned previously */
-
- v->sequence++;
+ /* Shift out any PCM/multipliers that we returned previously */
+ /* centerW is currently the center of the last block added */
if(v->pcm_returned && v->centerW>vi->blocksize[1]/2){
/* don't shift too much; we need to have a minimum PCM buffer of
1/2 long block */
- int shift=v->centerW-vi->blocksize[1]/2;
- shift=(v->pcm_returned<shift?v->pcm_returned:shift);
+ int shiftPCM=v->centerW-vi->blocksize[1]/2;
+ shiftPCM=(v->pcm_returned<shiftPCM?v->pcm_returned:shiftPCM);
- v->pcm_current-=shift;
- v->centerW-=shift;
- v->pcm_returned-=shift;
+ v->pcm_current-=shiftPCM;
+ v->centerW-=shiftPCM;
+ v->pcm_returned-=shiftPCM;
- if(shift){
+ if(shiftPCM){
int i;
for(i=0;i<vi->channels;i++)
- memmove(v->pcm[i],v->pcm[i]+shift,
+ memmove(v->pcm[i],v->pcm[i]+shiftPCM,
v->pcm_current*sizeof(double));
}
}
v->lW=v->W;
v->W=vb->W;
+ v->nW=-1;
v->gluebits+=vb->gluebits;
v->time_envelope_bits+=vb->time_envelope_bits;
int beginSl;
int endSl;
int i,j;
- double *windowL;
- double *windowN;
- /* Do we have enough PCM storage for the block? */
+ /* Do we have enough PCM/mult storage for the block? */
if(endW>v->pcm_storage){
- /* expand the PCM storage */
-
+ /* expand the storage */
v->pcm_storage=endW+vi->blocksize[1];
for(i=0;i<vi->channels;i++)
v->pcm[i]=realloc(v->pcm[i],v->pcm_storage*sizeof(double));
}
- /* Overlap/add */
+ /* overlap/add PCM */
+
switch(v->W){
case 0:
beginSl=0;
break;
}
- windowN=v->window[v->W][v->lW][v->lW];
- windowL=windowN+vi->blocksize[v->W]/2;
-
for(j=0;j<vi->channels;j++){
double *pcm=v->pcm[j]+beginW;
-
+
/* the overlap/add section */
for(i=beginSl;i<endSl;i++)
- pcm[i]=pcm[i]*windowL[i]+vb->pcm[j][i]*windowN[i];
+ pcm[i]+=vb->pcm[j][i];
/* the remaining section */
for(;i<sizeW;i++)
pcm[i]=vb->pcm[j][i];
function: codec headers
author: Monty <xiphmont@mit.edu>
modifications by: Monty
- last modification date: Oct 12 1999
+ last modification date: Oct 22 1999
********************************************************************/
typedef struct {
int winlen;
double *window;
+ mdct_lookup mdct;
} envelope_lookup;
typedef struct lpclook{
/* encode lookups */
- int *bscale;
+ int *uscale;
double *escale;
drft_lookup fft;
/* en/decode lookups */
int *iscale;
+ double *ifrac;
double *norm;
int n;
int ln;
int floororder[2];
int flooroctaves[2];
- int envelopesa;
- int envelopech;
int floorch;
/* int smallblock;
/* Encode-side settings for analysis */
+ int envelopesa;
double preecho_thresh;
double preecho_clamp;
logical bitstream */
int b_o_s; /* set after we've written the initial page
of a logical bitstream */
- long serialno;
- long pageno;
+ long serialno;
+ long pageno;
+ long packetno; /* sequence number for decode; the framing
+ knows where there's a hole in the data,
+ but we need coupling so that the codec
+ (which is in a seperate abstraction
+ layer) also knows about the gap */
int64_t pcmpos;
} ogg_stream_state;
long e_o_s;
int64_t frameno;
+ long packetno; /* sequence number for decode; the framing
+ knows where there's a hole in the data,
+ but we need coupling so that the codec
+ (which is in a seperate abstraction
+ layer) also knows about the gap */
} ogg_packet;
int pcm_current;
int pcm_returned;
- double **multipliers;
+ double *multipliers;
int envelope_storage;
int envelope_current;
typedef struct vorbis_block{
double **pcm;
- double **mult;
double **lpc;
double **lsp;
double *amp;
int pcm_channels; /* allocated, not used */
int pcm_storage; /* allocated, not used */
- int mult_channels; /* allocated, not used */
- int mult_storage; /* allocated, not used */
int floor_channels;
int floor_storage;
long W;
long nW;
int pcmend;
- int multend;
int eofflag;
int frameno;
+++ /dev/null
-#include <math.h>
-#include <stdio.h>
-#include "codec.h"
-
-#define READ 1024
-
-int main(){
- vorbis_dsp_state encode,decode;
- vorbis_info vi;
- vorbis_block vb;
- long counterin=0;
- long counterout=0;
- int done=0;
- char *temp[]={ "Test" ,"the Test band", "test records",NULL };
- int frame=0;
- int i;
-
- signed char buffer[READ*4+44];
-
- vi.channels=2;
- vi.rate=44100;
- vi.version=0;
- vi.mode=0;
- vi.user_comments=temp;
- vi.vendor="Xiphophorus";
- vi.smallblock=512;
- vi.largeblock=2048;
- vi.envelopesa=64;
- vi.envelopech=2;
-
- vi.floororder=30;
- vi.flooroctaves=5;
- vi.floorch=2;
-
- vi.balanceorder=4;
- vi.balanceoctaves=5;
-
- vi.channelmap=NULL;
- vi.preecho_thresh=10.;
- vi.preecho_clamp=.5;
-
- vorbis_analysis_init(&encode,&vi);
- vorbis_synthesis_init(&decode,&vi);
- vorbis_block_init(&encode,&vb);
-
- fread(buffer,1,44,stdin);
- fwrite(buffer,1,44,stdout);
-
- for(i=0;i<0;i++)
- fread(buffer,1,READ*4,stdin);
-
- while(!done){
- long bread=fread(buffer,1,READ*4,stdin);
- double **buf=vorbis_analysis_buffer(&encode,READ);
-
- /* uninterleave samples */
-
- for(i=0;i<bread/4;i++){
- buf[0][i]=((buffer[i*4+1]<<8)|(0x00ff&(int)buffer[i*4]))/32768.;
- buf[1][i]=((buffer[i*4+3]<<8)|(0x00ff&(int)buffer[i*4+2]))/32768.;
- }
-
- vorbis_analysis_wrote(&encode,i);
-
- while(vorbis_analysis_blockout(&encode,&vb)){
- double **pcm;
- int avail;
-
- /* analysis */
-
- vorbis_analysis(&vb);
- vorbis_analysis_packetout(&vb,&op);
-
- /* synthesis */
-
- vorbis_synthesis(&vb);
-
- counterin+=bread/4;
- vorbis_synthesis_blockin(&decode,&vb);
-
- while((avail=vorbis_synthesis_pcmout(&decode,&pcm))){
- if(avail>READ)avail=READ;
-
- for(i=0;i<avail;i++){
- int l=rint(pcm[0][i]*32768.);
- int r=rint(pcm[1][i]*32768.);
- if(l>32767)l=32767;
- if(r>32767)r=32767;
- if(l<-32768)l=-32768;
- if(r<-32768)r=-32768;
- buffer[i*4]=l&0xff;
- buffer[i*4+1]=(l>>8)&0xff;
- buffer[i*4+2]=r&0xff;
- buffer[i*4+3]=(r>>8)&0xff;
- }
-
- fwrite(buffer,1,avail*4,stdout);
-
- /*{
- FILE *out;
- char path[80];
- int i;
-
- sprintf(path,"syn%d",frame);
- out=fopen(path,"w");
-
- for(i=0;i<avail;i++)
- fprintf(out,"%ld %g\n",i+counterout,pcm[0][i]);
- fprintf(out,"\n");
- for(i=0;i<avail;i++)
- fprintf(out,"%ld %g\n",i+counterout,pcm[1][i]);
-
- fclose(out);
- }*/
-
- vorbis_synthesis_read(&decode,avail);
-
- counterout+=avail;
- frame++;
- }
-
-
- if(vb.eofflag){
- done=1;
- break;
- }
- }
- }
- return 0;
-}
function: PCM data envelope analysis and manipulation
author: Monty <xiphmont@mit.edu>
modifications by: Monty
- last modification date: Oct 17 1999
+ last modification date: Oct 22 1999
- Vorbis manipulates the dynamic range of the incoming PCM data
- envelope to minimise time-domain energy leakage from percussive and
- plosive waveforms being quantized in the MDCT domain.
+ Preecho calculation.
********************************************************************/
#include "mdct.h"
#include "envelope.h"
#include "bitwise.h"
+#include "window.h"
void _ve_envelope_init(envelope_lookup *e,int samples_per){
int i;
- e->winlen=samples_per*2;
+ e->winlen=samples_per;
e->window=malloc(e->winlen*sizeof(double));
+ mdct_init(&e->mdct,e->winlen);
+
+ /* We just use a straight sin(x) window for this */
+ for(i=0;i<e->winlen;i++)
+ e->window[i]=sin((i+.5)/e->winlen*M_PI);
- /* We just use a straight sin^2(x) window for this */
- for(i=0;i<e->winlen;i++){
- double temp=sin((i+.5)/e->winlen*M_PI);
- e->window[i]=temp*temp;
- }
}
void _ve_envelope_clear(envelope_lookup *e){
if(e->window)free(e->window);
+ mdct_clear(&e->mdct);
memset(e,0,sizeof(envelope_lookup));
}
-/* initial and final blocks are special cases. Eg:
- ______
- `--_
- |_______|_`-.___|_______|_______|
-
- ___
- _--' `--_
- |___.-'_|_______|_`-.___|_______|
-
- ___
- _--' `--_
- |_______|___.-'_|_______|_`-.___|
-
- _____
- _--'
- |_______|_______|____.-'|_______|
-
- as we go block by block, we watch the collective metrics span. If we
- span the threshhold (assuming the threshhold is active), we use an
- abbreviated vector */
-
-static int _ve_envelope_generate(double *mult,double *env,double *look,
- int n,int step){
- int i,j,p;
- double m,mo;
-
- for(j=0;j<n;j++)if(mult[j]!=1)break;
- if(j==n)return(0);
-
- n*=step;
- /* first multiplier special case */
- m=ldexp(1,mult[0]-1);
- for(p=0;p<step/2;p++)env[p]=m;
-
- /* mid multipliers normal case */
- for(j=1;p<n-step/2;j++){
- mo=m;
- m=ldexp(1,mult[j]-1);
- if(mo==m)
- for(i=0;i<step;i++,p++)env[p]=m;
- else
- for(i=0;i<step;i++,p++)env[p]=m*look[i]+mo*look[i+step];
- }
-
- /* last multiplier special case */
- for(;p<n;p++)env[p]=m;
- return(1);
-}
-
/* use MDCT for spectral power estimation */
-static void _ve_deltas(double *deltas,double *pcm,int n,double *win,
- int winsize){
+static void _ve_deltas(double *deltas,double *pcm,int n,double *window,
+ int winsize,mdct_lookup *m){
int i,j;
- mdct_lookup m;
double out[winsize/2];
-
- mdct_init(&m,winsize);
-
+
for(j=0;j<n;j++){
double acc=0.;
- mdct_forward(&m,pcm+j*winsize/2,out,win);
- for(i=0;i<winsize/2;i++)
+ mdct_forward(m,pcm+j*winsize,out,window);
+ for(i=winsize/10;i<winsize/2;i++)
acc+=fabs(out[i]);
if(deltas[j]<acc)deltas[j]=acc;
}
-
- mdct_clear(&m);
-
-#ifdef ANALYSIS
- {
- static int frameno=0;
- FILE *out;
- char buffer[80];
-
- sprintf(buffer,"delta%d.m",frameno++);
- out=fopen(buffer,"w+");
- for(j=0;j<n;j++)
- fprintf(out,"%g\n",deltas[j]);
- fclose(out);
- }
-#endif
}
-void _ve_envelope_multipliers(vorbis_dsp_state *v){
+void _ve_envelope_deltas(vorbis_dsp_state *v){
vorbis_info *vi=v->vi;
int step=vi->envelopesa;
- static int frame=0;
- /* we need a 1-1/4 envelope window overlap begin and 1/4 end */
- int dtotal=(v->pcm_current-step/2)/vi->envelopesa;
+ int dtotal=v->pcm_current/vi->envelopesa;
int dcurr=v->envelope_current;
- double *window=v->ve.window;
- int winlen=v->ve.winlen;
- int pch,ech;
+ int pch;
if(dtotal>dcurr){
- for(ech=0;ech<vi->envelopech;ech++){
- double *mult=v->multipliers[ech]+dcurr;
- memset(mult,0,sizeof(double)*(dtotal-dcurr));
+ double *mult=v->multipliers+dcurr;
+ memset(mult,0,sizeof(double)*(dtotal-dcurr));
- for(pch=0;pch<vi->channels;pch++){
-
- /* does this channel contribute to the envelope analysis */
- /*if(vi->envelopemap[pch]==ech){ not mapping yet */
- if(pch==ech){
-
- /* we need a 1/4 envelope window overlap front and back */
- double *pcm=v->pcm[pch]+dcurr*step-step/2;
- _ve_deltas(mult,pcm,dtotal-dcurr,window,winlen);
-
- }
- }
+ for(pch=0;pch<vi->channels;pch++){
+ double *pcm=v->pcm[pch]+dcurr*step;
+ _ve_deltas(mult,pcm,dtotal-dcurr,v->ve.window,v->ve.winlen,&v->ve.mdct);
}
v->envelope_current=dtotal;
- frame++;
- }
-}
-
-/* This readies the multiplier vector for use/coding. Clamp/adjust
- the multipliers to the allowed range and eliminate unneeded
- coefficients */
-
-void _ve_envelope_sparsify(vorbis_block *vb){
- vorbis_info *vi=vb->vd->vi;
- int ch;
- for(ch=0;ch<vi->envelopech;ch++){
- int flag=0;
- double *mult=vb->mult[ch];
- int n=vb->multend;
- double first=mult[0];
- double last=first;
- double clamp;
- int i;
#ifdef ANALYSIS
{
- static int frameno=0.;
+ static int frameno=0;
FILE *out;
- int j;
char buffer[80];
+ int j,k;
- sprintf(buffer,"del%d.m",frameno++);
+ sprintf(buffer,"delt%d.m",frameno);
out=fopen(buffer,"w+");
- for(j=0;j<n;j++)
- fprintf(out,"%g\n",mult[j]);
+ for(j=0;j<v->envelope_current;j++)
+ fprintf(out,"%d %g\n",j*vi->envelopesa+vi->envelopesa/2,v->multipliers[j]);
fclose(out);
- }
-#endif
- /* are we going to multiply anything? */
-
- for(i=1;i<n;i++){
- if(mult[i]>=last*vi->preecho_thresh){
- flag=1;
- break;
- }
- if(i<n-1 && mult[i+1]>=last*vi->preecho_thresh){
- flag=1;
- break;
- }
- last=mult[i];
- }
-
- if(flag){
- /* we need to adjust, so we might as well go nuts */
-
- int begin=i;
- clamp=last?last:1;
-
- for(i=0;i<begin;i++)mult[i]=0;
-
- for(;i<n;i++){
- if(mult[i]>=last*vi->preecho_thresh){
- last=mult[i];
-
- mult[i]=floor(log(mult[i]/clamp)/log(2));
- if(mult[i]>15)mult[i]=15;
- if(mult[i]<1)mult[i]=1;
-
- }else{
- mult[i]=0;
- }
- }
- }else
- memset(mult,0,sizeof(double)*n);
-
-#ifdef ANALYSIS
- {
- static int frameno=0.;
- FILE *out;
- int j;
- char buffer[80];
-
- sprintf(buffer,"sparse%d.m",frameno++);
+ sprintf(buffer,"deltpcm%d.m",frameno++);
out=fopen(buffer,"w+");
- for(j=0;j<n;j++)
- fprintf(out,"%g\n",mult[j]);
- fclose(out);
- }
-#endif
-
-
- }
-}
-
-void _ve_envelope_apply(vorbis_block *vb,int multp){
- static int frameno=0;
- vorbis_info *vi=vb->vd->vi;
- double env[vb->multend*vi->envelopesa];
- envelope_lookup *look=&vb->vd->ve;
- int i,j,k;
-
- for(i=0;i<vi->envelopech;i++){
- double *mult=vb->mult[i];
- double last=1.;
-
- /* fill in the multiplier placeholders */
-
- for(j=0;j<vb->multend;j++){
- if(mult[j]){
- last=mult[j];
- }else{
- mult[j]=last;
- }
- }
-
- /* compute the envelope curve */
- if(_ve_envelope_generate(mult,env,look->window,vb->multend,
- vi->envelopesa)){
-#ifdef ANALYSIS
- {
- FILE *out;
- char buffer[80];
-
- sprintf(buffer,"env%d.m",frameno);
- out=fopen(buffer,"w+");
- for(j=0;j<vb->multend*vi->envelopesa;j++)
- fprintf(out,"%g\n",env[j]);
- fclose(out);
- sprintf(buffer,"prepcm%d.m",frameno);
- out=fopen(buffer,"w+");
- for(j=0;j<vb->multend*vi->envelopesa;j++)
- fprintf(out,"%g\n",vb->pcm[i][j]);
- fclose(out);
- }
-#endif
-
for(k=0;k<vi->channels;k++){
-
- if(multp)
- for(j=0;j<vb->multend*vi->envelopesa;j++)
- vb->pcm[k][j]*=env[j];
- else
- for(j=0;j<vb->multend*vi->envelopesa;j++)
- vb->pcm[k][j]/=env[j];
-
- }
-
-#ifdef ANALYSIS
- {
- FILE *out;
- char buffer[80];
-
- sprintf(buffer,"pcm%d.m",frameno);
- out=fopen(buffer,"w+");
- for(j=0;j<vb->multend*vi->envelopesa;j++)
- fprintf(out,"%g\n",vb->pcm[i][j]);
- fclose(out);
+ for(j=0;j<v->pcm_current;j++)
+ fprintf(out,"%d %g\n",j,v->pcm[k][j]);
+ fprintf(out,"\n");
}
-#endif
+ fclose(out);
}
- frameno++;
- }
-}
-
-int _ve_envelope_encode(vorbis_block *vb){
- /* Not huffcoded yet. */
-
- vorbis_info *vi=vb->vd->vi;
- int scale=vb->W;
- int n=vb->multend;
- int i,j;
-
- /* range is 0-15 */
+#endif
- for(i=0;i<vi->envelopech;i++){
- double *mult=vb->mult[i];
- for(j=0;j<n;j++){
- _oggpack_write(&vb->opb,(int)(mult[j]),4);
- }
}
- return(0);
}
-/* synthesis expands the buffers in vb if needed. We can assume we
- have enough storage handed to us. */
-int _ve_envelope_decode(vorbis_block *vb){
- vorbis_info *vi=vb->vd->vi;
- int scale=vb->W;
- int n=vb->multend;
- int i,j;
-
- /* range is 0-15 */
-
- for(i=0;i<vi->envelopech;i++){
- double *mult=vb->mult[i];
- for(j=0;j<n;j++)
- mult[j]=_oggpack_read(&vb->opb,4);
- }
- return(0);
-}
function: PCM data envelope analysis and manipulation
author: Monty <xiphmont@mit.edu>
modifications by: Monty
- last modification date: Oct 07 1999
+ last modification date: Oct 22 1999
********************************************************************/
#ifndef _V_ENVELOPE_
#define _V_ENVELOPE_
-extern void _ve_envelope_multipliers(vorbis_dsp_state *v);
-extern void _ve_envelope_apply(vorbis_block *vb,int multp);
-extern void _ve_envelope_sparsify(vorbis_block *vb);
extern void _ve_envelope_init(envelope_lookup *e,int samples_per);
-extern void _ve_envelope_clear(envelope_lookup *e);
+extern void _ve_envelope_deltas(vorbis_dsp_state *v);
+extern void _ve_envelope_clear(envelope_lookup *e);
-extern int _ve_envelope_encode(vorbis_block *vb);
-extern int _ve_envelope_decode(vorbis_block *vb);
#endif
os->lacing_fill+=lacing_vals;
+ /* for the sake of completeness */
+ os->packetno++;
+
if(op->e_o_s)os->e_o_s=1;
return(0);
os->e_o_s=0;
os->b_o_s=0;
os->pageno=0;
+ os->packetno=0;
os->pcmpos=0;
return(0);
int ptr=os->lacing_returned;
if(os->lacing_packet<=ptr)return(0);
+
if(os->lacing_vals[ptr]&0x400){
/* We lost sync here; let the app know */
os->lacing_returned++;
+
+ /* we need to tell the codec there's a gap; it might need to
+ handle previous packet dependencies. */
+ os->packetno++;
return(-1);
}
if(val&0x200)op->e_o_s=0x200;
bytes+=size;
}
+
+ op->packetno=os->packetno;
op->frameno=os->pcm_vals[ptr];
op->bytes=bytes;
os->body_returned+=bytes;
os->lacing_returned=ptr+1;
}
+ os->packetno++;
return(1);
}
void checkpacket(ogg_packet *op,int len, int no, int pos){
long j;
+ static int sequence=0;
+ static int lastno=0;
+
if(op->bytes!=len){
fprintf(stderr,"incorrect packet length!\n");
exit(1);
exit(1);
}
+ /* packet number just follows sequence/gap; adjust the input number
+ for that */
+ if(no==0){
+ sequence=0;
+ }else{
+ sequence++;
+ if(no>lastno+1)
+ sequence++;
+ }
+ lastno=no;
+ if(op->packetno!=sequence){
+ fprintf(stderr,"incorrect packet sequence %ld != %d\n",op->packetno,sequence);
+ exit(1);
+ }
+
/* Test data */
for(j=0;j<op->bytes;j++)
if(op->packet[j]!=((j+no)&0xff)){
function: maintain the info structure, info <-> header packets
author: Monty <xiphmont@mit.edu>
modifications by: Monty
- last modification date: Oct 06 1999
+ last modification date: Oct 22 1999
********************************************************************/
memcpy(vi,&(predef_modes[mode]),sizeof(vorbis_info));
vi->threshhold_points=threshhold_points;
vi->user_comments=calloc(1,sizeof(char *));
- vi->vendor=strdup("Xiphophorus libVorbis I 19991018");
+ vi->vendor=strdup("Xiphophorus libVorbis I 19991022");
return(0);
}
vi->floororder[1]=_oggpack_read(&opb,8);
vi->flooroctaves[0]=_oggpack_read(&opb,8);
vi->flooroctaves[1]=_oggpack_read(&opb,8);
-
- vi->envelopesa=_oggpack_read(&opb,32);
- vi->envelopech=_oggpack_read(&opb,16);
vi->floorch=_oggpack_read(&opb,16);
return(0);
floor order for large block (octet)
floor octaves for small block (octet)
floor octaves for large block (octet)
-
- envelopesa (4 octets, lsb first)
- envelopech (4 octets, lsb first)
floorch (4 octets, lsb first)
*/
_oggpack_write(&opb,vi->floororder[1],8);
_oggpack_write(&opb,vi->flooroctaves[0],8);
_oggpack_write(&opb,vi->flooroctaves[1],8);
- _oggpack_write(&opb,vi->envelopesa,32);
- _oggpack_write(&opb,vi->envelopech,16);
_oggpack_write(&opb,vi->floorch,16);
/* build the packet */
function: LPC low level routines
author: Monty <monty@xiph.org>
modifications by: Monty
- last modification date: Oct 11 1999
+ last modification date: Oct 18 1999
********************************************************************/
#include "lpc.h"
#include "xlogmap.h"
-/* This is pared down for Vorbis where we only use LPC to encode
- spectral envelope curves. Thus we only are interested in
- generating the coefficients and recovering the curve from the
- coefficients. Autocorrelation LPC coeff generation algorithm
- invented by N. Levinson in 1947, modified by J. Durbin in 1959. */
+/* This is pared down for Vorbis. Autocorrelation LPC coeff generation
+ algorithm invented by N. Levinson in 1947, modified by J. Durbin in
+ 1959. */
-/* Input : n element envelope curve
+/* Input : n elements of time doamin data
Output: m lpc coefficients, excitation energy */
-double vorbis_gen_lpc(double *curve,double *lpc,lpc_lookup *l){
- int n=l->ln;
- int m=l->m;
- double aut[m+1],work[n+n],error;
- 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);
+double vorbis_lpc_from_data(double *data,double *lpc,int n,int m){
+ double aut[m+1],error;
+ int i,j;
- /* 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;
- }
-
/* autocorrelation, p+1 lag coefficients */
j=m+1;
while(j--){
double d=0;
- for(i=j;i<n;i++)d+=work[i]*work[i-j];
+ for(i=j;i<n;i++)d+=data[i]*data[i-j];
aut[j]=d;
}
-
+
/* Generate lpc coefficients from autocorr values */
error=aut[0];
r/=error;
/* Update LPC coefficients and total error */
-
+
lpc[i]=r;
for(j=0;j<i/2;j++){
double tmp=lpc[j];
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[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));
+}
+
/* On top of this basic LPC infrastructure we introduce two modifications:
1) Filter generation is limited in the resolution of features it
l->ln=mapped;
l->m=m;
l->iscale=malloc(n*sizeof(int));
+ l->ifrac=malloc(n*sizeof(double));
l->norm=malloc(n*sizeof(double));
for(i=0;i<n;i++){
if(encode_p){
/* encode */
- l->bscale=malloc(n*sizeof(int));
- l->escale=malloc(n*sizeof(double));
-
+ l->escale=malloc(mapped*sizeof(double));
+ l->uscale=malloc(n*sizeof(int));
+
+ /* undersample guard */
for(i=0;i<n;i++){
+ l->uscale[i]=rint(LOG_X(i,bias)/oct*mapped);
+ }
+
+ for(i=0;i<mapped;i++){
l->escale[i]=LINEAR_X(i/scale,bias);
- l->bscale[i]=rint(LOG_X(i,bias)*scale);
+ l->uscale[(int)(floor(l->escale[i]))]=-1;
+ l->uscale[(int)(ceil(l->escale[i]))]=-1;
}
+
}
/* decode; encode may use this too */
drft_init(&l->fft,mapped*2);
for(i=0;i<n;i++){
- l->iscale[i]=rint(LOG_X(i,bias)/oct*mapped);
- if(l->iscale[i]>=l->ln)l->iscale[i]=l->ln-1;
+ double is=LOG_X(i,bias)/oct*mapped;
+ if(is<0.)is=0.;
+
+ l->iscale[i]=floor(is);
+ if(l->iscale[i]>=l->ln-1)l->iscale[i]=l->ln-2;
+
+ l->ifrac[i]=is-floor(is);
+ if(l->ifrac[i]>1.)l->ifrac[i]=1.;
+
}
}
void lpc_clear(lpc_lookup *l){
if(l){
- if(l->bscale)free(l->bscale);
if(l->escale)free(l->escale);
drft_clear(&l->fft);
free(l->iscale);
+ free(l->ifrac);
free(l->norm);
}
}
}
- return vorbis_gen_lpc(work,lpc,l);
+ /* for(i=0;i<l->n;i++)
+ if(l->uscale[i]>0)
+ if(work[l->uscale[i]]<curve[i])work[l->uscale[i]]=curve[i];*/
+
+#ifdef ANALYSIS
+ {
+ int j;
+ FILE *out;
+ char buffer[80];
+ static int frameno=0;
+
+ 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);
}
int i;
_vlpc_de_helper(lcurve,lpc,amp,l);
+
+#ifdef ANALYSIS
+ {
+ int j;
+ FILE *out;
+ char buffer[80];
+ static int frameno=0;
+
+ sprintf(buffer,"loglpc%d.m",frameno++);
+ out=fopen(buffer,"w+");
+ for(j=0;j<l->ln;j++)
+ fprintf(out,"%g\n",lcurve[j]);
+ fclose(out);
+ }
+#endif
+
if(amp==0)return;
- for(i=0;i<l->n;i++)
- curve[i]=lcurve[l->iscale[i]]*l->norm[i];
+ for(i=0;i<l->n;i++){
+ int ii=l->iscale[i];
+ curve[i]=((1.-l->ifrac[i])*lcurve[ii]+
+ l->ifrac[i]*lcurve[ii+1])*l->norm[i];
+ }
+
}
void vorbis_lpc_apply(double *residue,double *lpc,double amp,lpc_lookup *l){
_vlpc_de_helper(lcurve,lpc,amp,l);
- for(i=0;i<l->n;i++)
- residue[i]*=lcurve[l->iscale[i]]*l->norm[i];
+ for(i=0;i<l->n;i++){
+ if(residue[i]!=0){
+ int ii=l->iscale[i];
+ residue[i]*=((1.-l->ifrac[i])*lcurve[ii]+
+ l->ifrac[i]*lcurve[ii+1])*l->norm[i];
+ }
+ }
+ }
+}
+
+/* subtract or add an lpc filter to data */
+
+void vorbis_lpc_residue(double *coeff,double *prime,int m,
+ double *data,long n){
+
+ /* in: coeff[0...m-1] LPC coefficients
+ prime[0...m-1] initial values
+ data[0...n-1] data samples
+ out: data[0...n-1] residuals from LPC prediction */
+
+ long i,j;
+ double work[m+n],y;
+
+ if(!prime)
+ for(i=0;i<m;i++)
+ work[i]=0;
+ else
+ for(i=0;i<m;i++)
+ work[i]=prime[i];
+
+ for(i=0;i<n;i++){
+ y=0;
+ for(j=0;j<m;j++)
+ y-=work[i+j]*coeff[m-j-1];
+
+ work[i+m]=data[i];
+ data[i]-=y;
+ }
+}
+
+
+void vorbis_lpc_predict(double *coeff,double *prime,int m,
+ double *data,long n){
+
+ /* in: coeff[0...m-1] LPC coefficients
+ prime[0...m-1] initial values (allocated size of n+m-1)
+ data[0...n-1] residuals from LPC prediction
+ out: data[0...n-1] data samples */
+
+ long i,j,o,p;
+ double y;
+ double work[n+m];
+
+ if(!prime)
+ for(i=0;i<m;i++)
+ work[i]=0.;
+ else
+ for(i=0;i<m;i++)
+ work[i]=prime[i];
+
+ for(i=0;i<n;i++){
+ y=data[i];
+ o=i;
+ p=m;
+ for(j=0;j<m;j++)
+ y-=work[o++]*coeff[--p];
+
+ data[i]=work[o]=y;
}
}
extern void lpc_clear(lpc_lookup *l);
/* simple linear scale LPC code */
-extern double vorbis_gen_lpc(double *curve,double *lpc,lpc_lookup *l);
-extern double vorbis_lpc_magnitude(double w,double *lpc, int m);
+extern double vorbis_lpc_from_data(double *data,double *lpc,int n,int m);
+extern double vorbis_lpc_from_spectrum(double *curve,double *lpc,lpc_lookup *l);
/* log scale layer */
extern double vorbis_curve_to_lpc(double *curve,double *lpc,lpc_lookup *l);
extern void vorbis_lpc_apply(double *residue,double *lpc, double amp,
lpc_lookup *l);
+/* standard lpc stuff */
+extern void vorbis_lpc_residue(double *coeff,double *prime,int m,
+ double *data,long n);
+extern void vorbis_lpc_predict(double *coeff,double *prime,int m,
+ double *data,long n);
+
#endif
}
}
-void mdct_backward(mdct_lookup *init, double *in, double *out){
+void mdct_backward(mdct_lookup *init, double *in, double *out, double *w){
int n=init->n;
double *x=alloca(n*sizeof(double));
int n2=n>>1;
int o3=n4+n2,o4=o3-1;
for(i=0;i<n4;i++){
- out[o1]=-(out[o2]=(*x * *BO - *(x+2) * *BE));
- out[o3]=out[o4]= -(*x * *BE + *(x+2) * *BO);
-
+ double temp1= (*x * *BO - *(x+2) * *BE);
+ double temp2=-(*x * *BE + *(x+2) * *BO);
+
+ out[o1]=-temp1*w[o1];
+ out[o2]= temp1*w[o2];
+ out[o3]= temp2*w[o3];
+ out[o4]= temp2*w[o4];
o1++;
o2--;
extern void mdct_clear(mdct_lookup *l);
extern void mdct_forward(mdct_lookup *init, double *in,
double *out, double *window);
-extern void mdct_backward(mdct_lookup *init, double *in, double *out);
+extern void mdct_backward(mdct_lookup *init, double *in,
+ double *out, double *window);
#endif
function: predefined encoding modes
author: Monty <xiphmont@mit.edu>
modifications by: Monty
- last modification date: Oct 15 1999
+ last modification date: Oct 22 1999
********************************************************************/
{ 2, 44100, 0, NULL, 0, NULL,
/* smallblock, largeblock, LPC order (small, large) */
{512, 2048}, {16,16},
- /* spectral octaves (small, large), envelope segment size */
- {5,5}, 64,
- /* envelope channels, spectrum channels */
- 1,2,
- /* preecho clamp trigger threshhold, clamp range, dummy */
- 4, 2, NULL,
+ /* spectral octaves (small, large), spectral channels */
+ {5,5}, 2,
+ /* thresh sample period, preecho clamp trigger threshhold, range, dummy */
+ 128, 4, 2, NULL,
/* noise masking curve dB attenuation levels [20] */
- {-12,-12,-18,-18,-18,-18,-18,-10,-5,-2,
- 0,0,0,0,1,2,3,3,4,5},
+ {-12,-12,-18,-18,-18,-18,-18,-18,-18,-12,
+ -8,-4,0,0,1,2,3,3,4,5},
+ /*{-100,-100,-100,-100,-100,-100,-100,-24,-24,-24,
+ -24,-24,-24,-24,-24,-24,-24,-24,-24,-24}*/
/* noise masking scale biases */
- .95,1.01,.001,
+ .95,1.01,.01,
/* tone masking curve dB attenuation levels [20] */
- {-24,-24,-24,-24,-24,-24,-24,-24,-24,-24,
- -24,-24,-24,-24,-24,-24,-24,-24,-24,-24},
+ {-20,-20,-20,-20,-20,-20,-20,-20,-20,-20,
+ -20,-20,-20,-20,-20,-20,-20,-20,-20,-20},
/* tone masking rolloff settings (dB per octave), octave bias */
90,60,.001,
NULL,NULL,NULL},
}
}
-/*************************************************************************/
-/* Continuous balance analysis/synthesis *********************************/
-
-
-/* Compute the average continual spectral balance of the given vectors
- (in radians); encode into LPC coefficients */
-
-double _vp_balance_compute(double *A, double *B, double *lpc,lpc_lookup *vb){
- /* correlate in time (the response function is small). Log
- frequency scale, small mapping */
-
- int n=vb->n;
- int mapped=vb->ln;
-
- /* 256/15 are arbitrary but unimportant to decoding */
- int resp=15;
- int i;
-
- /* This is encode side. Don't think too hard about it */
-
- double workA[mapped+resp];
- double workB[mapped+resp];
- double p[mapped];
- double workC[resp];
-
- memset(workA,0,sizeof(workA));
- memset(workB,0,sizeof(workB));
- memset(workC,0,sizeof(workC));
-
- for(i=0;i<n;i++){
- int j=vb->bscale[i]+(resp>>1);
- double mag_sq=A[i]*A[i]+B[i]*B[i];
- double phi;
-
- if(B[i]==0)
- phi=M_PI/2.;
- else{
- phi=atan(A[i]/B[i]);
- if((A[i]<0 && B[i]>0)||
- (A[i]>0 && B[i]<0)){
- /* rotate II and IV into the first quadrant. III is already there */
- phi+=M_PI/2;
- }
- }
-
- workA[j]+=mag_sq*sin(phi);
- workB[j]+=mag_sq*cos(phi);
- }
-
- /* prepare convolution vector. Play with a few different shapes */
-
- for(i=0;i<resp;i++){
- workC[i]=sin(M_PI*(i+1)/(float)(resp+1));
- workC[i]*=workC[i];
- }
-
- time_convolve(workA,workC,mapped,resp);
- time_convolve(workB,workC,mapped,resp);
-
- {
- double amp1;
-
- for(i=0;i<mapped;i++){
- p[i]=atan2(workA[i],workB[i]);
- }
-
- amp1=sqrt(vorbis_gen_lpc(p,lpc,vb));
-
- return(amp1);
- }
-
-}
-
-/*void _vp_balance_apply(double *A, double *B, double *lpc, double amp,
- lpc_lookup *vb,int divp){
- int i;
- for(i=0;i<vb->n;i++){
- double mag=sqrt(A[i]*A[i]+B[i]*B[i]);
- double del=vorbis_lpc_magnitude(vb->dscale[i],lpc,vb->m)*amp;
- double phi=atan2(A[i],B[i]);
-
- if(divp)
- phi-=del;
- else
- phi+=del;
-
- A[i]=mag*sin(phi);
- B[i]=mag*cos(phi);
- }
- }*/
int scale=vb->W;
int m=vb->vd->vi->floororder[scale];
- int n=vb->pcmend/2;
+ int n=vb->pcmend*4;
int last=0;
double dlast=0.;
double min=M_PI/n/2.;
int _vs_spectrum_decode(vorbis_block *vb,double *amp,double *lsp){
int scale=vb->W;
int m=vb->vd->vi->floororder[scale];
- int n=vb->pcmend/2;
+ int n=vb->pcmend*4;
int last=0;
double dlast=0.;
int bits=rint(log(n)/log(2));
for(i=0;i<m;i++){
int val=_oggpack_read(&vb->opb,bits);
lsp[i]=(last+=val)*M_PI/n;
+
/* Underpowered but sufficient */
if(lsp[i]<dlast+min)lsp[i]=dlast+min;
dlast=lsp[i];
if(val>16)val=16;
if(val<-16)val=-16;
- if(val==0 || val==2 || val==-2){
+ /*if(val==0 || val==2 || val==-2){
if(data[i]<0){
val=-1;
}else{
val=1;
}
- }
+ }*/
data[i]=val;
/*if(val<0){
#include "spectrum.h"
int vorbis_synthesis(vorbis_block *vb,ogg_packet *op){
+ double *window;
vorbis_dsp_state *vd=vb->vd;
vorbis_info *vi=vd->vi;
oggpack_buffer *opb=&vb->opb;
return(-1);
}
- /* Encode the block size */
+ /* Decode the block size */
vb->W=_oggpack_read(opb,1);
+ if(vb->W){
+ vb->lW=_oggpack_read(opb,1);
+ vb->nW=_oggpack_read(opb,1);
+ }else{
+ vb->lW=0;
+ vb->nW=0;
+ }
+
+ window=vb->vd->window[vb->W][vb->lW][vb->nW];
/* other random setup */
vb->frameno=op->frameno;
+ vb->sequence=op->packetno-3; /* first block is third packet */
+
vb->eofflag=op->e_o_s;
vl=&vb->vd->vl[vb->W];
spectral_order=vi->floororder[vb->W];
/* The storage vectors are large enough; set the use markers */
n=vb->pcmend=vi->blocksize[vb->W];
- vb->multend=n/vi->envelopesa;
- /* recover the time envelope */
- if(_ve_envelope_decode(vb)<0)return(-1);
+ /* No envelope encoding yet */
+ _oggpack_write(opb,0,1);
for(i=0;i<vi->channels;i++){
double *lpc=vb->lpc[i];
/* recover the spectral residue */
if(_vs_residue_decode(vb,vb->pcm[i])<0)return(-1);
+#ifdef ANALYSIS
+ {
+ int j;
+ FILE *out;
+ char buffer[80];
+
+ sprintf(buffer,"Sres%d.m",vb->sequence);
+ out=fopen(buffer,"w+");
+ for(j=0;j<n/2;j++)
+ fprintf(out,"%g\n",vb->pcm[i][j]);
+ fclose(out);
+ }
+#endif
+
/* LSP->LPC */
vorbis_lsp_to_lpc(lsp,lpc,vl->m);
/* apply envelope to residue */
+#ifdef ANALYSIS
+ {
+ int j;
+ FILE *out;
+ char buffer[80];
+ double curve[n/2];
+ vorbis_lpc_to_curve(curve,lpc,vb->amp[i],vl);
+
+
+ sprintf(buffer,"Smask%d.m",vb->sequence);
+ out=fopen(buffer,"w+");
+ for(j=0;j<n/2;j++)
+ fprintf(out,"%g\n",curve[j]);
+ fclose(out);
+
+ sprintf(buffer,"Slsp%d.m",vb->sequence);
+ out=fopen(buffer,"w+");
+ for(j=0;j<vl->m;j++)
+ fprintf(out,"%g\n",lsp[j]);
+ fclose(out);
+
+ sprintf(buffer,"Slpc%d.m",vb->sequence);
+ out=fopen(buffer,"w+");
+ for(j=0;j<vl->m;j++)
+ fprintf(out,"%g\n",lpc[j]);
+ fclose(out);
+ }
+#endif
+
vorbis_lpc_apply(vb->pcm[i],lpc,vb->amp[i],vl);
+
+#ifdef ANALYSIS
+ {
+ int j;
+ FILE *out;
+ char buffer[80];
+
+ sprintf(buffer,"Sspectrum%d.m",vb->sequence);
+ out=fopen(buffer,"w+");
+ for(j=0;j<n/2;j++)
+ fprintf(out,"%g\n",vb->pcm[i][j]);
+ fclose(out);
+ }
+#endif
+
/* MDCT->time */
- mdct_backward(&vb->vd->vm[vb->W],vb->pcm[i],vb->pcm[i]);
+ mdct_backward(&vb->vd->vm[vb->W],vb->pcm[i],vb->pcm[i],window);
}
-
- /* apply time domain envelope */
- _ve_envelope_apply(vb,1);
-
return(0);
}
projects / platform / upstream / libvorbis.git / commitdiff