lib/analysis.c

   1 /********************************************************************
   2  *                                                                  *
   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.                            *
   7  *                                                                  *
   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/                                             *
  11  *                                                                  *
  12  ********************************************************************
  13
  14  function: single-block PCM analysis
  15  last mod: $Id: analysis.c,v 1.20 2000/01/05 03:10:53 xiphmont Exp $
  16
  17  ********************************************************************/
  18
  19 #include <stdio.h>
  20 #include <string.h>
  21 #include <math.h>
  22 #include "vorbis/codec.h"
  23
  24 #include "os.h"
  25 #include "lpc.h"
  26 #include "lsp.h"
  27 #include "envelope.h"
  28 #include "mdct.h"
  29 #include "psy.h"
  30 #include "bitwise.h"
  31 #include "spectrum.h"
  32
  33 /* this code is still seriously abbreviated.  I'm filling in pieces as
  34    we go... --Monty 19991004 */
  35
  36 int vorbis_analysis(vorbis_block *vb,ogg_packet *op){
  37   int i;
  38   double           *window=vb->vd->window[vb->W][vb->lW][vb->nW];
  39   psy_lookup       *vp=&vb->vd->vp[vb->W];
  40   lpc_lookup       *vl=&vb->vd->vl[vb->W];
  41   vorbis_dsp_state *vd=vb->vd;
  42   vorbis_info      *vi=vd->vi;
  43   oggpack_buffer   *opb=&vb->opb;
  44
  45   int              n=vb->pcmend;
  46   int              spectral_order=vi->floororder[vb->W];
  47
  48   vb->gluebits=0;
  49   vb->time_envelope_bits=0;
  50   vb->spectral_envelope_bits=0;
  51   vb->spectral_residue_bits=0;
  52
  53   /*lpc_lookup       *vbal=&vb->vd->vbal[vb->W];
  54     double balance_v[vbal->m];
  55     double balance_amp;*/
  56
  57   /* first things first.  Make sure encode is ready*/
  58   _oggpack_reset(opb);
  59   /* Encode the packet type */
  60   _oggpack_write(opb,0,1);
  61
  62   /* Encode the block size */
  63   _oggpack_write(opb,vb->W,1);
  64   if(vb->W){
  65     _oggpack_write(opb,vb->lW,1);
  66     _oggpack_write(opb,vb->nW,1);
  67   }
  68
  69   /* No envelope encoding yet */
  70   _oggpack_write(opb,0,1);
  71
  72   /* time domain PCM -> MDCT domain */
  73   for(i=0;i<vi->channels;i++)
  74     mdct_forward(&vd->vm[vb->W],vb->pcm[i],vb->pcm[i],window);
  75
  76   /* no balance yet */
  77
  78   /* extract the spectral envelope and residue */
  79   /* just do by channel.  No coupling yet */
  80   {
  81     for(i=0;i<vi->channels;i++){
  82       static int frameno=0;
  83       int j;
  84       double *floor=alloca(n/2*sizeof(double));
  85       double *curve=alloca(n/2*sizeof(double));
  86       double *lpc=vb->lpc[i];
  87       double *lsp=vb->lsp[i];
  88
  89       memset(floor,0,sizeof(double)*n/2);
  90
  91 #ifdef ANALYSIS
  92       {
  93         FILE *out;
  94         char buffer[80];
  95
  96         sprintf(buffer,"Aspectrum%d.m",vb->sequence);
  97         out=fopen(buffer,"w+");
  98         for(j=0;j<n/2;j++)
  99           fprintf(out,"%g\n",vb->pcm[i][j]);
 100         fclose(out);
 101
 102       }
 103 #endif
 104
 105       _vp_mask_floor(vp,vb->pcm[i],floor);
 106
 107 #ifdef ANALYSIS
 108       {
 109         FILE *out;
 110         char buffer[80];
 111
 112         sprintf(buffer,"Apremask%d.m",vb->sequence);
 113         out=fopen(buffer,"w+");
 114         for(j=0;j<n/2;j++)
 115           fprintf(out,"%g\n",floor[j]);
 116         fclose(out);
 117       }
 118 #endif
 119
 120       /* Convert our floor to a set of lpc coefficients */
 121       vb->amp[i]=sqrt(vorbis_curve_to_lpc(floor,lpc,vl));
 122
 123       /* LSP <-> LPC is orthogonal and LSP quantizes more stably */
 124       vorbis_lpc_to_lsp(lpc,lsp,vl->m);
 125
 126       /* code the spectral envelope; mutates the lsp coeffs to reflect
 127          what was actually encoded */
 128       _vs_spectrum_encode(vb,vb->amp[i],lsp);
 129
 130       /* Generate residue from the decoded envelope, which will be
 131          slightly different to the pre-encoding floor due to
 132          quantization.  Slow, yes, but perhaps more accurate */
 133
 134       vorbis_lsp_to_lpc(lsp,lpc,vl->m);
 135       vorbis_lpc_to_curve(curve,lpc,vb->amp[i],vl);
 136
 137       /* this may do various interesting massaging too...*/
 138       if(vb->amp[i])_vs_residue_train(vb,vb->pcm[i],curve,n/2);
 139       _vs_residue_quantize(vb->pcm[i],curve,vi,n/2);
 140
 141 #ifdef ANALYSIS
 142       {
 143         FILE *out;
 144         char buffer[80];
 145
 146         sprintf(buffer,"Alpc%d.m",vb->sequence);
 147         out=fopen(buffer,"w+");
 148         for(j=0;j<vl->m;j++)
 149           fprintf(out,"%g\n",lpc[j]);
 150         fclose(out);
 151
 152         sprintf(buffer,"Alsp%d.m",vb->sequence);
 153         out=fopen(buffer,"w+");
 154         for(j=0;j<vl->m;j++)
 155           fprintf(out,"%g\n",lsp[j]);
 156         fclose(out);
 157
 158         sprintf(buffer,"Amask%d.m",vb->sequence);
 159         out=fopen(buffer,"w+");
 160         for(j=0;j<n/2;j++)
 161           fprintf(out,"%g\n",curve[j]);
 162         fclose(out);
 163
 164         sprintf(buffer,"Ares%d.m",vb->sequence);
 165         out=fopen(buffer,"w+");
 166         for(j=0;j<n/2;j++)
 167           fprintf(out,"%g\n",vb->pcm[i][j]);
 168         fclose(out);
 169       }
 170 #endif
 171
 172       /* encode the residue */
 173       _vs_residue_encode(vb,vb->pcm[i]);
 174
 175     }
 176   }
 177
 178   /* set up the packet wrapper */
 179
 180   op->packet=opb->buffer;
 181   op->bytes=_oggpack_bytes(opb);
 182   op->b_o_s=0;
 183   op->e_o_s=vb->eofflag;
 184   op->frameno=vb->frameno;
 185   op->packetno=vb->sequence; /* for sake of completeness */
 186
 187   return(0);
 188 }
 189
 190
 191
 192
 193 /* commented out, relocated balance stuff */
 194   /*{
 195     double *C=vb->pcm[0];
 196     double *D=vb->pcm[1];
 197
 198     balance_amp=_vp_balance_compute(D,C,balance_v,vbal);
 199
 200     {
 201       FILE *out;
 202       char buffer[80];
 203
 204       sprintf(buffer,"com%d.m",frameno);
 205       out=fopen(buffer,"w+");
 206       for(i=0;i<n/2;i++){
 207         fprintf(out," 0. 0.\n");
 208         fprintf(out,"%g %g\n",C[i],D[i]);
 209         fprintf(out,"\n");
 210       }
 211       fclose(out);
 212
 213       sprintf(buffer,"L%d.m",frameno);
 214       out=fopen(buffer,"w+");
 215       for(i=0;i<n/2;i++){
 216         fprintf(out,"%g\n",C[i]);
 217       }
 218       fclose(out);
 219       sprintf(buffer,"R%d.m",frameno);
 220       out=fopen(buffer,"w+");
 221       for(i=0;i<n/2;i++){
 222         fprintf(out,"%g\n",D[i]);
 223       }
 224       fclose(out);
 225
 226     }
 227
 228     _vp_balance_apply(D,C,balance_v,balance_amp,vbal,1);
 229
 230     {
 231       FILE *out;
 232       char buffer[80];
 233
 234       sprintf(buffer,"bal%d.m",frameno);
 235       out=fopen(buffer,"w+");
 236       for(i=0;i<n/2;i++){
 237         fprintf(out," 0. 0.\n");
 238         fprintf(out,"%g %g\n",C[i],D[i]);
 239         fprintf(out,"\n");
 240       }
 241       fclose(out);
 242       sprintf(buffer,"C%d.m",frameno);
 243       out=fopen(buffer,"w+");
 244       for(i=0;i<n/2;i++){
 245         fprintf(out,"%g\n",C[i]);
 246       }
 247       fclose(out);
 248       sprintf(buffer,"D%d.m",frameno);
 249       out=fopen(buffer,"w+");
 250       for(i=0;i<n/2;i++){
 251         fprintf(out,"%g\n",D[i]);
 252       }
 253       fclose(out);
 254
 255     }
 256   }*/