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: Oct 21 1999
  18
  19  ********************************************************************/
  20
  21 #include <stdio.h>
  22 #include <string.h>
  23 #include <math.h>
  24 #include "lpc.h"
  25 #include "lsp.h"
  26 #include "codec.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[n/2];
  85       double curve[n/2];
  86       double *lpc=vb->lpc[i];
  87       double *lsp=vb->lsp[i];
  88
  89       memset(floor,0,sizeof(double)*n/2);
  90
  91       _vp_noise_floor(vp,vb->pcm[i],floor);
  92
  93 #ifdef ANALYSIS
  94       {
  95         FILE *out;
  96         char buffer[80];
  97
  98         sprintf(buffer,"Aspectrum%d.m",vb->sequence);
  99         out=fopen(buffer,"w+");
 100         for(j=0;j<n/2;j++)
 101           fprintf(out,"%g\n",vb->pcm[i][j]);
 102         fclose(out);
 103
 104         sprintf(buffer,"Anoise%d.m",vb->sequence);
 105         out=fopen(buffer,"w+");
 106         for(j=0;j<n/2;j++)
 107           fprintf(out,"%g\n",floor[j]);
 108         fclose(out);
 109       }
 110 #endif
 111
 112       _vp_mask_floor(vp,vb->pcm[i],floor);
 113
 114 #ifdef ANALYSIS
 115       {
 116         FILE *out;
 117         char buffer[80];
 118
 119         sprintf(buffer,"Apremask%d.m",vb->sequence);
 120         out=fopen(buffer,"w+");
 121         for(j=0;j<n/2;j++)
 122           fprintf(out,"%g\n",floor[j]);
 123         fclose(out);
 124       }
 125 #endif
 126
 127       /* Convert our floor to a set of lpc coefficients */
 128       vb->amp[i]=sqrt(vorbis_curve_to_lpc(floor,lpc,vl));
 129
 130       /* LSP <-> LPC is orthogonal and LSP quantizes more stably */
 131       vorbis_lpc_to_lsp(lpc,lsp,vl->m);
 132
 133       /* code the spectral envelope; mutates the lsp coeffs to reflect
 134          what was actually encoded */
 135       _vs_spectrum_encode(vb,vb->amp[i],lsp);
 136
 137       /* Generate residue from the decoded envelope, which will be
 138          slightly different to the pre-encoding floor due to
 139          quantization.  Slow, yes, but perhaps more accurate */
 140
 141       vorbis_lsp_to_lpc(lsp,lpc,vl->m);
 142       vorbis_lpc_to_curve(curve,lpc,vb->amp[i],vl);
 143
 144       /* this may do various interesting massaging too...*/
 145       _vs_residue_quantize(vb->pcm[i],curve,vi,n/2);
 146
 147 #ifdef ANALYSIS
 148       {
 149         FILE *out;
 150         char buffer[80];
 151
 152         sprintf(buffer,"Alpc%d.m",vb->sequence);
 153         out=fopen(buffer,"w+");
 154         for(j=0;j<vl->m;j++)
 155           fprintf(out,"%g\n",lpc[j]);
 156         fclose(out);
 157
 158         sprintf(buffer,"Alsp%d.m",vb->sequence);
 159         out=fopen(buffer,"w+");
 160         for(j=0;j<vl->m;j++)
 161           fprintf(out,"%g\n",lsp[j]);
 162         fclose(out);
 163
 164         sprintf(buffer,"Amask%d.m",vb->sequence);
 165         out=fopen(buffer,"w+");
 166         for(j=0;j<n/2;j++)
 167           fprintf(out,"%g\n",curve[j]);
 168         fclose(out);
 169
 170         sprintf(buffer,"Ares%d.m",vb->sequence);
 171         out=fopen(buffer,"w+");
 172         for(j=0;j<n/2;j++)
 173           fprintf(out,"%g\n",vb->pcm[i][j]);
 174         fclose(out);
 175       }
 176 #endif
 177
 178       /* encode the residue */
 179       _vs_residue_encode(vb,vb->pcm[i]);
 180
 181     }
 182   }
 183
 184   /* set up the packet wrapper */
 185
 186   op->packet=opb->buffer;
 187   op->bytes=_oggpack_bytes(opb);
 188   op->b_o_s=0;
 189   op->e_o_s=vb->eofflag;
 190   op->frameno=vb->frameno;
 191   op->packetno=vb->sequence; /* for sake of completeness */
 192
 193   return(0);
 194 }
 195
 196
 197
 198
 199 /* commented out, relocated balance stuff */
 200   /*{
 201     double *C=vb->pcm[0];
 202     double *D=vb->pcm[1];
 203
 204     balance_amp=_vp_balance_compute(D,C,balance_v,vbal);
 205
 206     {
 207       FILE *out;
 208       char buffer[80];
 209
 210       sprintf(buffer,"com%d.m",frameno);
 211       out=fopen(buffer,"w+");
 212       for(i=0;i<n/2;i++){
 213         fprintf(out," 0. 0.\n");
 214         fprintf(out,"%g %g\n",C[i],D[i]);
 215         fprintf(out,"\n");
 216       }
 217       fclose(out);
 218
 219       sprintf(buffer,"L%d.m",frameno);
 220       out=fopen(buffer,"w+");
 221       for(i=0;i<n/2;i++){
 222         fprintf(out,"%g\n",C[i]);
 223       }
 224       fclose(out);
 225       sprintf(buffer,"R%d.m",frameno);
 226       out=fopen(buffer,"w+");
 227       for(i=0;i<n/2;i++){
 228         fprintf(out,"%g\n",D[i]);
 229       }
 230       fclose(out);
 231
 232     }
 233
 234     _vp_balance_apply(D,C,balance_v,balance_amp,vbal,1);
 235
 236     {
 237       FILE *out;
 238       char buffer[80];
 239
 240       sprintf(buffer,"bal%d.m",frameno);
 241       out=fopen(buffer,"w+");
 242       for(i=0;i<n/2;i++){
 243         fprintf(out," 0. 0.\n");
 244         fprintf(out,"%g %g\n",C[i],D[i]);
 245         fprintf(out,"\n");
 246       }
 247       fclose(out);
 248       sprintf(buffer,"C%d.m",frameno);
 249       out=fopen(buffer,"w+");
 250       for(i=0;i<n/2;i++){
 251         fprintf(out,"%g\n",C[i]);
 252       }
 253       fclose(out);
 254       sprintf(buffer,"D%d.m",frameno);
 255       out=fopen(buffer,"w+");
 256       for(i=0;i<n/2;i++){
 257         fprintf(out,"%g\n",D[i]);
 258       }
 259       fclose(out);
 260
 261     }
 262   }*/