hsigp.c

该压缩包为最新版htk的源代码,htk是现在比较流行的语音处理软件,请有兴趣的朋友下载使用

/* ----------------------------------------------------------- */
/*                                                             */
/*                          ___                                */
/*                       |_| | |_/   SPEECH                    */
/*                       | | | | \   RECOGNITION               */
/*                       =========   SOFTWARE                  */ 
/*                                                             */
/*                                                             */
/* ----------------------------------------------------------- */
/* developed at:                                               */
/*                                                             */
/*      Speech Vision and Robotics group                       */
/*      Cambridge University Engineering Department            */
/*      http://svr-www.eng.cam.ac.uk/                          */
/*                                                             */
/*      Entropic Cambridge Research Laboratory                 */
/*      (now part of Microsoft)                                */
/*                                                             */
/* ----------------------------------------------------------- */
/*         Copyright: Microsoft Corporation                    */
/*          1995-2000 Redmond, Washington USA                  */
/*                    http://www.microsoft.com                 */
/*                                                             */
/*              2001  Cambridge University                     */
/*                    Engineering Department                   */
/*                                                             */
/*   Use of this software is governed by a License Agreement   */
/*    ** See the file License for the Conditions of Use  **    */
/*    **     This banner notice must not be removed      **    */
/*                                                             */
/* ----------------------------------------------------------- */
/*      File: HSigP.c:   Signal Processing Routines            */
/* ----------------------------------------------------------- */

char *hsigp_version = "!HVER!HSigP:   3.3 [CUED 28/04/05]";
char *hsigp_vc_id = "$Id: HSigP.c,v 1.1.1.1 2005/05/12 10:52:51 jal58 Exp $";

#include "HShell.h"        /* HTK Libraries */
#include "HMem.h"
#include "HMath.h"
#include "HSigP.h"

/*
   This module provides a set of basic speech signal processing
   routines and feature level transformations.
*/

/* ------------------------ Trace Flags ------------------------- */

static int trace = 0;
#define T_MEL  0002     /* Mel filterbank */

/* -------------------- Config and Memory ----------------------- */

static MemHeap sigpHeap;
static ConfParam *cParm[MAXGLOBS];       /* config parameters */
static int numParm = 0;

/* ---------------------- Initialisation -------------------------*/

/* EXPORT->InitSigP: initialise the SigP module */
void InitSigP(void)
{
   int i;

   Register(hsigp_version,hsigp_vc_id);
   numParm = GetConfig("HSIGP", TRUE, cParm, MAXGLOBS);
   if (numParm>0){
      if (GetConfInt(cParm,numParm,"TRACE",&i)) trace = i;
   }
   CreateHeap(&sigpHeap,"sigpHeap",MSTAK,1,0.0,5000,5000);
}

/* --------------- Windowing and PreEmphasis ---------------------*/

/* ZeroMean: zero mean a complete speech waveform */
void ZeroMean(short *data, long nSamples)
{
   long i,hiClip=0,loClip=0;
   short *x;
   double sum=0.0,y,mean;

   x = data;
   for (i=0;i<nSamples;i++,x++)
      sum += *x;
   mean = sum / (float)nSamples;
   x = data;
   for (i=0;i<nSamples;i++,x++){
      y = (double)(*x) - mean;
      if (y<-32767.0){
         y = -32767.0; ++loClip;
      }
      if (y>32767.0){
         y = 32767.0; ++hiClip;
      }
      *x = (short) ((y>0.0) ? y+0.5 : y-0.5);
   }
   if (loClip>0)
      HError(-5322,"ZeroMean: %d samples too -ve\n",loClip);
   if (hiClip>0)
      HError(-5322,"ZeroMean: %d samples too +ve\n",hiClip);
}

static int hamWinSize = 0;          /* Size of current Hamming window */
static Vector hamWin = NULL;        /* Current Hamming window */

/* GenHamWindow: generate precomputed Hamming window function */
static void GenHamWindow (int frameSize)
{
   int i;
   float a;
   
   if (hamWin==NULL || VectorSize(hamWin) < frameSize)
      hamWin = CreateVector(&sigpHeap,frameSize);
   a = TPI / (frameSize - 1);
   for (i=1;i<=frameSize;i++)
      hamWin[i] = 0.54 - 0.46 * cos(a*(i-1));
   hamWinSize = frameSize;
}

/* EXPORT->Ham: Apply Hamming Window to Speech frame s */
void Ham (Vector s)
{
   int i,frameSize;
   
   frameSize=VectorSize(s);
   if (hamWinSize != frameSize)
      GenHamWindow(frameSize);
   for (i=1;i<=frameSize;i++)
      s[i] *= hamWin[i];
}

/* EXPORT->PreEmphasise: pre-emphasise signal in s */
void PreEmphasise (Vector s, float k)
{
   int i;
   float preE;
   
   preE = k;
   for (i=VectorSize(s);i>=2;i--)
      s[i] -= s[i-1]*preE;
   s[1] *= 1.0-preE;
}

/* --------------- Linear Prediction Coding Operations ------------- */

/* AutoCorrelate: store auto coef 1 to p in r and return energy r[0] */
static float AutoCorrelate(Vector s, Vector r, int p, int frameSize)
{
   float sum,energy;
   int   i,j;

   energy = 0.0;
   for (i=0;i<=p;i++) {
      sum = 0.0;
      for (j=1;j<=frameSize-i;j++)
         sum += s[j]*s[j+i];
      if (i==0) 
         energy = sum;
      else
         r[i] = sum;
   }
   return energy;
}

/* Durbins recursion to get LP coeffs for auto values */
static float Durbin(Vector k, Vector thisA, Vector r, float E, int order)
{
   Vector newA;
   float ki;         /* Current Reflection Coefficient */
   int i,j;
 
   newA  = CreateVector(&gstack,order);
   for (i=1;i<=order;i++) {
      ki = r[i];              /* Calc next reflection coef */
      for (j=1;j<i;j++)
         ki = ki + thisA[j] * r[i - j];
      ki = ki / E;   
      if (k!=NULL) k[i] = ki;
      E *= 1 - ki*ki;         /* Update Error */
      newA[i] = -ki;          /* Calc new filter coef */
      for (j=1;j<i;j++)
         newA[j] = thisA[j] - ki * thisA[i - j];
      for (j=1;j<=i;j++)   
         thisA[j] = newA[j];
   }
   FreeVector(&gstack,newA);
   return (E);
}

/* EXPORT->Wave2LPC: Calculate LPCoef in a & RefC in k */
void Wave2LPC (Vector s, Vector a, Vector k, float *re, float *te)
{
   Vector thisA;     /* Current LP filter coefficients */
   Vector r;         /* AutoCorrelation Sequence */
   float E;          /* Prediction Error */
   int   p,frameSize;

   if (a==NULL && k==NULL)
      HError(5320,"Wave2LPC: Null a and k vectors in WaveToLPC");  
   if (a!=NULL) 
      p=VectorSize(a); 
   else
      p=VectorSize(k);
   r = CreateVector(&gstack,p);
   thisA = (a!=NULL)?a:CreateVector(&gstack,p);
   frameSize=VectorSize(s);
   E = AutoCorrelate(s,r,p,frameSize);
   *te = E;
   *re = Durbin(k,thisA,r,E,p);
   FreeVector(&gstack,r);
}

/* EXPORT->LPC2RefC: transfer from filter to ref coef */
void LPC2RefC(Vector a, Vector k)
{
   Vector thisA; /* Current LP filter coefficients */
   Vector newA;  /* New LP filter coefficients */
   int i,j,p;
   float ki,x;
   
   p=VectorSize(a);
   thisA = CreateVector(&gstack,p);
   newA  = CreateVector(&gstack,p);
   CopyVector(a,thisA);
   for (i=p;i>=1;i--)  { 
      ki = -thisA[i];
      k[i] = ki;
      x = 1 - ki*ki;
      for (j=1;j<i;j++) 
         newA[j] = (thisA[j] + ki * thisA[i - j]) / x;
      for (j=1;j<i;j++) 
         thisA[j] = newA[j];
   }
   FreeVector(&gstack,thisA);
}

/* EXPORT->RefC2LPC: transfer from ref coef to filter */
void RefC2LPC (Vector k, Vector a)
{
   Vector thisA; /* Current LP filter coefficients */
   Vector newA;  /* New LP filter coefficients */
   int i,j,p;
   float ki;
   
   p=VectorSize(k);
   thisA = CreateVector(&gstack,p);
   newA  = CreateVector(&gstack,p);
   for (i=1;i<=p;i++) { 
      ki = k[i];
      newA[i] = -ki;
      for (j=1;j<i;j++) 
         newA[j] = thisA[j] - ki * thisA[i - j];
      for (j=1;j<=i;j++)   
         thisA[j] = newA[j];
   }
   for (i=1;i<=p;i++)  a[i]=thisA[i];
   FreeVector(&gstack,thisA);
}

/* EXPORT->LPC2Cepstrum: transfer from lpc to cepstral coef */
void LPC2Cepstrum (Vector a, Vector c)
{
   int i,n,p;
   float sum;
   
   p=VectorSize(c);
   for (n=1;n<=p;n++)  { 
      sum = 0.0;
      for (i=1;i<n;i++) 
         sum = sum + (n - i) * a[i] * c[n - i];
      c[n] = -(a[n] + sum / n);
   }
}

/* EXPORT->Cepstrum2LPC: transfer from cepstral coef to lpc */
void Cepstrum2LPC (Vector c, Vector a)
{
   int i,n,p;
   float sum;
   
   p=VectorSize(a);
   for (n=1;n<=p;n++)  { 
      sum = 0.0;
      for (i=1;i<n;i++) 
         sum = sum + (n - i) * a[i] * c[n - i];
      a[n] = -(c[n] + sum / n);
   }
}

/* -------------------- FFT Based Operations ----------------------- */

/* EXPORT-> FVec2Spectrum: cvt feature vector f to a spectrum, fzero
   is the value of the 0'th feature vector coefficient which
   is typically omitted by HSigP routines eg a0 = 1.0 for LPC
*/
void FVec2Spectrum (float fzero, Vector f, Vector s)
{
   int i,p,n;
   
   p=VectorSize(f); n=VectorSize(s);
   s[1] = fzero;
   for (i=1;i<=p;i++) 
      s[i+1] = f[i];
   for (i=p+2;i<=n;i++) 
      s[i] = 0.0;
   Realft(s);
}

/* EXPORT-> FFT: apply fft/invfft to complex s */
void FFT(Vector s, int invert)
{
   int ii,jj,n,nn,limit,m,j,inc,i;
   double wx,wr,wpr,wpi,wi,theta;
   double xre,xri,x;
   
   n=VectorSize(s);
   nn=n / 2; j = 1;
   for (ii=1;ii<=nn;ii++) {
      i = 2 * ii - 1;
      if (j>i) {
         xre = s[j]; xri = s[j + 1];
         s[j] = s[i];  s[j + 1] = s[i + 1];
         s[i] = xre; s[i + 1] = xri;
      }
      m = n / 2;
      while (m >= 2  && j > m) {
         j -= m; m /= 2;
      }
      j += m;
   };
   limit = 2;
   while (limit < n) {
      inc = 2 * limit; theta = TPI / limit;
      if (invert) theta = -theta;
      x = sin(0.5 * theta);
      wpr = -2.0 * x * x; wpi = sin(theta); 
      wr = 1.0; wi = 0.0;
      for (ii=1; ii<=limit/2; ii++) {
         m = 2 * ii - 1;
         for (jj = 0; jj<=(n - m) / inc;jj++) {
            i = m + jj * inc;
            j = i + limit;
            xre = wr * s[j] - wi * s[j + 1];
            xri = wr * s[j + 1] + wi * s[j];
            s[j] = s[i] - xre; s[j + 1] = s[i + 1] - xri;
            s[i] = s[i] + xre; s[i + 1] = s[i + 1] + xri;
         }
         wx = wr;
         wr = wr * wpr - wi * wpi + wr;
         wi = wi * wpr + wx * wpi + wi;
      }
      limit = inc;
   }
   if (invert)
      for (i = 1;i<=n;i++) 
         s[i] = s[i] / nn;
   
}

/* EXPORT-> Realft: apply fft to real s */
void Realft (Vector s)
{
   int n, n2, i, i1, i2, i3, i4;
   double xr1, xi1, xr2, xi2, wrs, wis;
   double yr, yi, yr2, yi2, yr0, theta, x;

   n=VectorSize(s) / 2; n2 = n/2;
   theta = PI / n;
   FFT(s, FALSE);
   x = sin(0.5 * theta);
   yr2 = -2.0 * x * x;
   yi2 = sin(theta); yr = 1.0 + yr2; yi = yi2;
   for (i=2; i<=n2; i++) {
      i1 = i + i - 1;      i2 = i1 + 1;
      i3 = n + n + 3 - i2; i4 = i3 + 1;
      wrs = yr; wis = yi;
      xr1 = (s[i1] + s[i3])/2.0; xi1 = (s[i2] - s[i4])/2.0;
      xr2 = (s[i2] + s[i4])/2.0; xi2 = (s[i3] - s[i1])/2.0;
      s[i1] = xr1 + wrs * xr2 - wis * xi2;
      s[i2] = xi1 + wrs * xi2 + wis * xr2;
      s[i3] = xr1 - wrs * xr2 + wis * xi2;
      s[i4] = -xi1 + wrs * xi2 + wis * xr2;
      yr0 = yr;
      yr = yr * yr2 - yi  * yi2 + yr;
      yi = yi * yr2 + yr0 * yi2 + yi;
   }
   xr1 = s[1];
   s[1] = xr1 + s[2];
   s[2] = 0.0;
}
   
/* EXPORT-> SpecModulus: store modulus of s in m */
void SpecModulus(Vector s, Vector m)
{
   int i,j;
   float x,y;
   
   for (i=1;i<=VectorSize(s)/2;i++) {
      j=i+i; x=s[j-1]; y=s[j];
      m[i]=sqrt(x*x + y*y);
   }
}

/* EXPORT-> SpecLogModulus: store log modulus of s in m */
void SpecLogModulus(Vector s, Vector m, Boolean invert)
{
   int i,j;
   float x,y;
   
   for (i=1;i<=VectorSize(s)/2;i++) {
      j=i+i; x=s[j-1]; y=s[j];
      x=0.5*log(x*x + y*y);
      m[i] = invert ? -x : x;
   }
}

/* EXPORT-> SpecPhase: store phase of s in m */
void SpecPhase(Vector s, Vector m)
{
   int i,j;
   float ph,re,im;
   
   for (i=1;i<=VectorSize(s)/2;i++) {
      j=i+i;
      re=s[j-1]; im=s[j];
      if (re==0.0) 
         ph = (im>=0.0) ? PI/2.0 : -PI/2.0;
      else {
         ph=atan(im/re);
         if (ph<0.0 && re<0.0)
            ph += PI;
         else if (ph>0.0 && im<0.0)
            ph -= PI;
      }
      m[i]=ph;
   }
}

/* -------------------- MFCC Related Operations -------------------- */

/* EXPORT->Mel: return mel-frequency corresponding to given FFT index */
float Mel(int k,float fres)
{
   return 1127 * log(1 + (k-1)*fres);
}

/* EXPORT->WarpFreq: return warped frequency */
float WarpFreq (float fcl, float fcu, float freq, float minFreq, float maxFreq , float alpha)
{
   if (alpha == 1.0)
      return freq;
   else {
      float scale = 1.0 / alpha;
      float cu = fcu * 2 / (1 + scale);
      float cl = fcl * 2 / (1 + scale);

      float au = (maxFreq - cu * scale) / (maxFreq - cu);
      float al = (cl * scale - minFreq) / (cl - minFreq);
      
      if (freq > cu)
         return  au * (freq - cu) + scale * cu ;
      else if (freq < cl)
         return al * (freq - minFreq) + minFreq ;
      else
         return scale * freq ;
   }