Spectrum.cpp

Go to the documentation of this file.

/**********************************************************************
 
  Audacity: A Digital Audio Editor
 
  Spectrum.cpp
 
  Dominic Mazzoni
 
*******************************************************************//*******************************************************************/
 
#include "Spectrum.h"
 
#include <math.h>
 
#include "SampleFormat.h"
 
bool ComputeSpectrum(const float * data, size_t width,
                     size_t windowSize,
                     double WXUNUSED(rate), float *output,
                     bool autocorrelation, int windowFunc)
{
   if (width < windowSize)
      return false;
 
   if (!data || !output)
      return true;
 
   Floats processed{ windowSize };
 
   for (size_t i = 0; i < windowSize; i++)
      processed[i] = float(0.0);
   auto half = windowSize / 2;
 
   Floats in{ windowSize };
   Floats out{ windowSize };
   Floats out2{ windowSize };
 
   size_t start = 0;
   unsigned windows = 0;
   while (start + windowSize <= width) {
      for (size_t i = 0; i < windowSize; i++)
         in[i] = data[start + i];
 
      WindowFunc(windowFunc, windowSize, in.get());
 
      if (autocorrelation) {
         // Take FFT
         RealFFT(windowSize, in.get(), out.get(), out2.get());
         // Compute power
         for (size_t i = 0; i < windowSize; i++)
            in[i] = (out[i] * out[i]) + (out2[i] * out2[i]);
 
         // Tolonen and Karjalainen recommend taking the cube root
         // of the power, instead of the square root
 
         for (size_t i = 0; i < windowSize; i++)
            in[i] = powf(in[i], 1.0f / 3.0f);
 
         // Take FFT
         RealFFT(windowSize, in.get(), out.get(), out2.get());
      }
      else
         PowerSpectrum(windowSize, in.get(), out.get());
 
      // Take real part of result
      for (size_t i = 0; i < half; i++)
        processed[i] += out[i];
 
      start += half;
      windows++;
   }
 
   if (autocorrelation) {
 
      // Peak Pruning as described by Tolonen and Karjalainen, 2000
      /*
       Combine most of the calculations in a single for loop.
       It should be safe, as indexes refer only to current and previous elements,
       that have already been clipped, etc...
      */
      for (size_t i = 0; i < half; i++) {
        // Clip at zero, copy to temp array
        if (processed[i] < 0.0)
            processed[i] = float(0.0);
        out[i] = processed[i];
        // Subtract a time-doubled signal (linearly interp.) from the original
        // (clipped) signal
        if ((i % 2) == 0)
           processed[i] -= out[i / 2];
        else
           processed[i] -= ((out[i / 2] + out[i / 2 + 1]) / 2);
 
        // Clip at zero again
        if (processed[i] < 0.0)
            processed[i] = float(0.0);
      }
 
      // Reverse and scale
      for (size_t i = 0; i < half; i++)
         in[i] = processed[i] / (windowSize / 4);
      for (size_t i = 0; i < half; i++)
         processed[half - 1 - i] = in[i];
   } else {
      // Convert to decibels
      // But do it safely; -Inf is nobody's friend
      for (size_t i = 0; i < half; i++){
         float temp=(processed[i] / windowSize / windows);
         if (temp > 0.0)
            processed[i] = 10 * log10(temp);
         else
            processed[i] = 0;
      }
   }
 
   for(size_t i = 0; i < half; i++)
      output[i] = processed[i];
 
   return true;
}