title: Cepstrum and Formant Tracking Using Apple Accelerate Framework
body: <p>I've been using this web page as a guideline for formant tracking of speech...</p>
  
  <p><a href="http://iitg.vlab.co.in/?sub=59&amp;brch=164&amp;sim=615&amp;cnt=1" rel="nofollow noreferrer">http://iitg.vlab.co.in/?sub=59&amp;brch=164&amp;sim=615&amp;cnt=1</a></p>
  
  <p>It all seems to be going pretty well, except for the last step, which is the converting of the cepstrum into a smoothed representation for simple peak picking for the formant tracking.  The spectrograph looks good, and the cepstrograph (can I say that? :P) also looks good (from what I can tell), but the final stage the results (smoothed formant representation) are not what I expected.</p>
  
  <p>I uploaded a sample of each stage as visual images to...</p>
  
  <p><a href="https://imgur.com/a/62duS" rel="nofollow noreferrer">http://imgur.com/a/62duS</a></p>
  
  <p>This sample is for the speech of the sound 'i' as in 'beed'.  According to this site...</p>
  
  <p><a href="http://home.cc.umanitoba.ca/~robh/howto.html#formants" rel="nofollow noreferrer">http://home.cc.umanitoba.ca/~robh/howto.html#formants</a></p>
  
  <p>the first formant should come in around 500hz, and the second and third around 2200hz and 2800 hz respectively.  The spetrograph shows something very similar, but on the last stage I am gettings results similar to...</p>
  
  <p>F1 - 891
  F2 - 1550
  F3 - 2329</p>
  
  <p>Any insight would be greatly appreciated.  I've been going round in circles on this for some time.  My code looks as follows...</p>
  
  <pre><code>// set up fft parameters
  UInt32 log2n = 9;
  UInt32 n = 512;
  UInt32 window = n;
  UInt32 halfN = n/2;
  UInt32 stride = 1;
  FFTSetup setupReal = [appDelegate getFftSetup];
  int stepSize = (hpBuffer.sampleCount-window) / quantizeCount;
  
  // calculate volume from raw samples, because it seems more reliable that fft
  UInt32 volumeWindow = 128;
  volumeBuffer = malloc(sizeof(float)*quantizeCount);
  int windowPos = 0;
  for (int i=0; i &lt; quantizeCount; i++) {
      windowPos += stepSize;
      float total = 0.0f;
      float max = 0.0f;
      for (int p=windowPos; p &lt; windowPos+volumeWindow; p++) {
          total += sampleBuffer.buffer[p];
          if (sampleBuffer.buffer[p] &gt; max)
              max = sampleBuffer.buffer[p];
      }
      volumeBuffer[i] = max;
  }
  
  // normalize volumebuffer
  [FloatAudioBuffer normalizePositiveBuffer:volumeBuffer ofSize:quantizeCount];
  
  // allocate memory for complex array
  COMPLEX_SPLIT complexArray;
  complexArray.realp = (float*)malloc(4096*sizeof(float));
  complexArray.imagp = (float*)malloc(4096*sizeof(float));
  
  // allocate some space for temporary hamming buffer
  float *hamBuffer = malloc(n*sizeof(float));
  
  // create spectrum and feature buffer
  spectrumBuffer = malloc(sizeof(float)*halfN*quantizeCount);
  formantBuffer = malloc(sizeof(float)*4096*quantizeCount);
  cepstrumBuffer = malloc(sizeof(float)*halfN*quantizeCount);
  lowCepstrumBuffer = malloc(sizeof(float)*featureCount*quantizeCount);
  featureBuffer = malloc(sizeof(float)*featureCount*quantizeCount);
  
  // create data point for each quantize segment
  float TWOPI = 2.0f * M_PI;
  for (int s=0; s &lt; quantizeCount; s++) {
      // copy buffer data into a seperate array and apply hamming window
      int offset = (int)(s * stepSize);
      for (int i=0; i &lt; n; i++)
          hamBuffer[i] = hpBuffer.buffer[offset+i] * ((1.0f-0.46f) - 0.46f*cos(TWOPI*i/((float)n-1.0f)));
  
      // configure float array into acceptable input array format (interleaved)
      vDSP_ctoz((COMPLEX*)hamBuffer, 2, &amp;complexArray, 1, halfN);
  
      // run FFT
      vDSP_fft_zrip(setupReal, &amp;complexArray, stride, log2n, FFT_FORWARD);
  
      // Absolute square (equivalent to mag^2)
      complexArray.imagp[0] = 0.0f;
      vDSP_zvmags(&amp;complexArray, 1, complexArray.realp, 1, halfN);
      bzero(complexArray.imagp, (halfN) * sizeof(float));
  
      // scale
      float scale = 1.0f / (2.0f*(float)n);
      vDSP_vsmul(complexArray.realp, 1, &amp;scale, complexArray.realp, 1, halfN);
  
      // get log of absolute values for passing to inverse FFT for cepstrum
      for (int i=0; i &lt; halfN; i++)
          complexArray.realp[i] = logf(sqrtf(complexArray.realp[i]));
  
      // save this into spectrum buffer
      memcpy(&amp;spectrumBuffer[s*halfN], complexArray.realp, halfN*sizeof(float));
  
  
      // convert spectrum to interleaved ready for inverse fft
      vDSP_ctoz((COMPLEX*)&amp;spectrumBuffer[s*halfN], 2, &amp;complexArray, 1, halfN/2);
  
      // create cepstrum
      vDSP_fft_zrip(setupReal, &amp;complexArray, stride, log2n-1, FFT_INVERSE);
  
      //convert interleaved to real and straight into cepstrum buffer
      vDSP_ztoc(&amp;complexArray, 1, (COMPLEX*)&amp;cepstrumBuffer[s*halfN], 2, halfN/2);
  
      // copy first part of cepstrum into low cepstrum buffer
      memcpy(&amp;lowCepstrumBuffer[s*featureCount], &amp;cepstrumBuffer[s*halfN], featureCount*sizeof(float));
  
      // make 8000 point array based on the first 15 values
      float *tempArray = malloc(8192*sizeof(float));
      for (int i=0; i &lt; 8192; i++) {
          if (i &lt; 15)
              tempArray[i] = cepstrumBuffer[s*halfN+i];
          else
              tempArray[i] = 0.0f;
      }
      vDSP_ctoz((COMPLEX*)tempArray, 2, &amp;complexArray, 1, 4096);
      float newLog2n = log2f(8192.0f);
      complexArray.imagp[0] = 0.0f;
      vDSP_fft_zrip(setupReal, &amp;complexArray, stride, newLog2n, FFT_FORWARD);
      vDSP_zvmags(&amp;complexArray, 1, complexArray.realp, 1, 4096);
      bzero(complexArray.imagp, (4096) * sizeof(float));
  
      // scale
      scale = 1.0f / (2.0f*(float)8192);
      vDSP_vsmul(complexArray.realp, 1, &amp;scale, complexArray.realp, 1, 4096);
  
      // get magnitude
      for (int i=0; i &lt; 4096; i++)
          complexArray.realp[i] = sqrtf(complexArray.realp[i]);
  
      // write to formant buffer
      memcpy(&amp;formantBuffer[s*4096], complexArray.realp, 4096*sizeof(float));
  
      // complex array now contains formant spectrum
      // it's large, so get features here!
      // try simple peak picking algorithm for first 3 formants
      int formantIndex = 0;
      float *peaks = malloc(6*sizeof(float));
      for (int i=0; i &lt; 6; i++)
          peaks[i] = 0.0f;
      for (int i=1; i &lt; 4096-1 &amp;&amp; formantIndex &lt; 6; i++) {
          if (complexArray.realp[i-1] &lt; complexArray.realp[i] &amp;&amp;
              complexArray.realp[i+1] &lt; complexArray.realp[i])
              peaks[formantIndex++] = i;
      }
  </code></pre>