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