view src/feature_extractor.cpp @ 0:add35537fdbb tip

Initial import

//  Copyright 2011, Ian Hobson.
//
//  This file is part of gpsynth.
//
//  gpsynth is free software: you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation, either version 3 of the License, or
//  (at your option) any later version.
//
//  gpsynth is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with gpsynth in the file COPYING. 
//  If not, see http://www.gnu.org/licenses/.

#include "feature_extractor.hpp"

#include "statistics.hpp"
#include "std_ex.hpp"
#include "window_makers.hpp"

#include "sndfile.hh"

#include "boost/bind.hpp"
#include "boost/filesystem.hpp"

#include <algorithm>
#include <cmath>
#include <iostream>
#include <stdexcept>
#include <sstream>

namespace bf = boost::filesystem;

namespace dsp {

FeatureExtractor::FeatureExtractor(const std::string& file_path /* = "" */,
                                   int window_size /* = 2048 */,
                                   int hop_size /* = 512 */) 
: window_size_(-1), // initialized by SetFrameSize
  hop_size_(hop_size),
  bin_frequencies_(window_size),
  spectrum_analyzer_(window_size),
  mfcc_analyzer_(window_size, 40, 13, 30, 10000)
{
  SetWindowSize(window_size);
  if (!file_path.empty()) {
    LoadFile(file_path);
  }
}

void FeatureExtractor::ClearCaches() {
  spectrum_.clear();
  magnitude_spectrum_.clear();
  log_magnitude_spectrum_.clear();
  spectral_centroid_.clear();
  spectral_spread_.clear();
  spectral_flux_.clear();
  pitch_.clear();
  energy_.clear();
  mfccs_.clear();
  delta_mfccs_.clear();
  double_delta_mfccs_.clear();
}
  
void FeatureExtractor::SetWindowSize(int window_size) {
  if (window_size_ != window_size) {
    window_size_ = window_size;
    frame_.resize(window_size);
    // generate new window
    window_.resize(window_size);
    dsp::HammingWindow(window_size, window_.begin());
    ClearCaches();
  }
}
  
void FeatureExtractor::SetHopSize(int hop_size) {
  if (hop_size_ != hop_size) {
    hop_size_ = hop_size;
    ClearCaches();
  }
}
  
void FeatureExtractor::LoadFile(const std::string& file_path) {
  SndfileHandle file(file_path);
  if (file.error()) {
    std::stringstream message;
    message << "FeatureExtractor::ExtractFeaturesFromFile: "
            << "Error reading from file '" << file_path << "' "
            << '(' << file.error() << ')';
    throw std::runtime_error(message.str());
  }
  sf_count_t frames = file.frames();
  int channels = file.channels();
  file_data_.resize(frames * channels);
  // read the audio data in to the buffer
  file.readf(&file_data_[0], frames);
  // only supporting mono at the moment, remove extra channels
  if (channels > 1) {
    for (int i = 0; i < frames; i++) {
      file_data_[i] = file_data_[i * channels];
    }
    file_data_.resize(frames);
  }
  // success, clear caches
  ClearCaches();
  // update file info
  file_path_ = file_path;
  sample_rate_ = file.samplerate();
  frames_ = static_cast<int>(file_data_.size()) / hop_size_;
  if (frames_ == 0) {
    frames_++;
  }
  // bin info
  bin_size_ = sample_rate_ / static_cast<Value>(window_size_);
  std::generate(bin_frequencies_.begin(), bin_frequencies_.end(),
                stdx::Stepper<Value>(bin_size_));
}
  
Value FeatureExtractor::Duration() const {
  return file_data_.size() / static_cast<Value>(sample_rate_);
}

const std::vector<ComplexValueList>& FeatureExtractor::Spectrum() {
  if (spectrum_.empty()) {
    const int frame_size = window_size_;
    const int half_frame_size = frame_size / 2;
    const int hop_size = hop_size_;
    spectrum_.resize(frames_);
    ValueList::const_iterator data_begin = file_data_.begin();
    ValueList::const_iterator data_end = file_data_.end();
    for (int i = 0; i < frames_; i++) {
      //std::cout << "frame:" << i << '\n';
      // find the frame bounds
      ValueList::const_iterator read_start = data_begin + (i * hop_size);
      ValueList::const_iterator read_end = read_start + frame_size;
      if (read_end > data_end) {
        read_end = data_end;
      }
      // extract the frame and apply window
      std::transform(read_start, read_end, window_.begin(), 
                     frame_.begin(), std::multiplies<Value>());
      // on the last frame we'll need to fill the rest of the frame with zeros
      std::size_t data_read = read_end - read_start;
      if (data_read < frame_size) {
        std::fill(frame_.begin() + data_read, frame_.end(), 0);
      }
      // get the spectrum of the current frame
      std::copy(frame_.begin(), frame_.end(), spectrum_analyzer_.Input());
      spectrum_analyzer_.Execute();
      spectrum_[i].assign(spectrum_analyzer_.Output(),
                          spectrum_analyzer_.Output() + half_frame_size);
    }
  }
  return spectrum_;
}

const std::vector<ValueList>& FeatureExtractor::MagnitudeSpectrum() {
  if (magnitude_spectrum_.empty()) {
    // get complex spectrum
    const std::vector<ComplexValueList>& spectrum = Spectrum();
    // get magnitude spectrum for each frame
    std::size_t frames = spectrum.size();
    magnitude_spectrum_.resize(frames);
    for (int i = 0; i < frames; i++) {
      // take the absolute value of each complex bin
      magnitude_spectrum_[i].resize(spectrum[i].size());
      std::transform(spectrum[i].begin(), spectrum[i].end(),
                     magnitude_spectrum_[i].begin(), std::abs<Value>);
    }
  }
  return magnitude_spectrum_;
}
  
const std::vector<ValueList>& FeatureExtractor::LogMagnitudeSpectrum() {
  if (log_magnitude_spectrum_.empty()) {
    // get magnitude spectrum
    const std::vector<ValueList>& magnitude_spectrum = MagnitudeSpectrum();
    // get magnitude spectrum for each frame
    std::size_t frames = magnitude_spectrum.size();
    log_magnitude_spectrum_.resize(frames);
    for (int i = 0; i < frames; i++) {
      log_magnitude_spectrum_[i].resize(magnitude_spectrum[i].size());
      // copy magnitudes to log magnitude buffer, avoiding zeros
      std::replace_copy(magnitude_spectrum[i].begin(),
                        magnitude_spectrum[i].end(), 
                        log_magnitude_spectrum_[i].begin(),
                        0.0,
                        std::numeric_limits<Value>::min());
      // take the log of each magnitude (ugly casting for cmath function...)
      std::transform(log_magnitude_spectrum_[i].begin(),
                     log_magnitude_spectrum_[i].end(),
                     log_magnitude_spectrum_[i].begin(),
                     static_cast<double(*)(double)>(std::log10));
    }
  }
  return log_magnitude_spectrum_;
}

const ValueList& FeatureExtractor::SpectralCentroid() {
  if (spectral_centroid_.empty()) {
    MagnitudeSpectrum();
    std::size_t frames = magnitude_spectrum_.size();
    spectral_centroid_.resize(frames);
    std::size_t bins = magnitude_spectrum_[0].size();
    ValueList weighted_frequencies(bins);
    for (int i = 0; i < frames; i++) {
      const ValueList& magnitudes = magnitude_spectrum_[i];
      // get bin frequencies weighted by bin magnitudes
      std::transform(magnitudes.begin(),
                     magnitudes.end(),
                     bin_frequencies_.begin(),
                     weighted_frequencies.begin(),
                     std::multiplies<Value>());
      // divide mean of weighted frequencies by mean of bin magnitudes
      Value magnitude = stats::Sum(magnitudes);
      if (magnitude > 0) {
        spectral_centroid_[i] = stats::Sum(weighted_frequencies) / magnitude;
      } else {
        spectral_centroid_[i] = 0;
      }
    }
  }
  return spectral_centroid_;
}


const ValueList& FeatureExtractor::SpectralSpread() {
  if (spectral_spread_.empty()) {
    SpectralCentroid();
    std::size_t frames = magnitude_spectrum_.size();
    spectral_spread_.resize(frames);
    std::size_t bins = magnitude_spectrum_[0].size();
    ValueList weighted_difference(bins);
    for (int i = 0; i < frames; i++) {
      Value centroid = spectral_centroid_[i];
      const ValueList& magnitudes = magnitude_spectrum_[i];
      for (int bin = 0; bin < bins; bin++) {
        Value difference = bin_frequencies_[bin] - centroid;
        weighted_difference[bin] = difference * difference * magnitudes[bin];
      }
      Value magnitude = stats::Sum(magnitudes);
      if (magnitude > 0) {
        Value variance = stats::Sum(weighted_difference) / magnitude;
        spectral_spread_[i] = std::sqrt(variance);
      } else {
        spectral_spread_[i] = 0;
      }
    }
  }
  return spectral_spread_;
}

const ValueList& FeatureExtractor::Pitch() {
  if (pitch_.empty()) {
    MagnitudeSpectrum();
    pitch_.resize(frames_);
    std::deque<Value> recent_pitches;
    const int moving_average_size = 8;
    for (int i = 0; i < frames_; i++) {
      const ValueList& mags = magnitude_spectrum_[i];
      // find the peak bin number 
      std::ptrdiff_t peak_bin = std::distance(mags.begin(),
                                              std::max_element(mags.begin(),
                                                               mags.end()));
      // interpolate peak frequency
      double offset;
      double centre = mags[peak_bin];
      if (centre == 0) {
        pitch_[i] = 0;
      } else {
        if (peak_bin == 0) {
          // parabolic interpolation of first 2 bin magnitudes
          double above = mags[peak_bin + 1];
          offset = above / (2 * (2 * centre - above));        
        } else if (peak_bin == mags.size() - 1) {
          // parabolic interpolation of last 2 bin magnitudes
          double below = mags[peak_bin - 1];
          offset = (-below) / (2 * (2 * centre - below));
        } else {
          double below = mags[peak_bin - 1];
          double above = mags[peak_bin + 1];
          // gaussian interpolation
          offset = std::log(above / below) 
                   / (2 * std::log((centre * centre) / (above * below)));
        }
        Value frame_pitch = bin_size_ * (peak_bin + offset);
        if (recent_pitches.size() == moving_average_size) {
          recent_pitches.pop_front();
        }
        recent_pitches.push_back(frame_pitch);
        pitch_[i] = stats::Mean(recent_pitches);
      }
    }
  }
  return pitch_;
}
  
const ValueList& FeatureExtractor::Energy() {
  if (energy_.empty()) {
    const int frame_size = window_size_;
    const int hop_size = hop_size_;
    energy_.resize(frames_);
    ValueList::const_iterator data_begin = file_data_.begin();
    ValueList::const_iterator data_end = file_data_.end();
    for (int i = 0; i < frames_; i++) {
      // find the frame bounds
      ValueList::const_iterator read_start = data_begin + (i * hop_size);
      ValueList::const_iterator read_end = read_start + frame_size;
      if (read_end > data_end) {
        read_end = data_end;
      }
      Value energy(0);
      std::ptrdiff_t samples = std::distance(read_start, read_end);
      while (read_start != read_end) {
        energy += std::pow(*read_start++, 2);
      }
      energy_[i] = energy / samples;
    }
  }
  return energy_;
}

const std::vector<ValueList>& FeatureExtractor::MFCCs() {
  if (mfccs_.empty()) {
    const std::vector<ValueList>& magnitude_spectrum = MagnitudeSpectrum();
    mfcc_analyzer_.SetSampleRate(sample_rate_);
    mfccs_.resize(frames_);
    for (int i = 0; i < frames_; i++) {
      mfccs_[i] = mfcc_analyzer_.ExtractMFCCS(magnitude_spectrum[i]);
    }
  }
  return mfccs_;
}
  
const std::vector<ValueList>& FeatureExtractor::DeltaMFCCs() {
  if (delta_mfccs_.empty()) {
    const std::vector<ValueList>& mfccs = MFCCs();
    // resize buffers
    const std::size_t number_of_mfccs = mfccs[0].size();
    ValueList empty_values(number_of_mfccs, 0);
    delta_mfccs_.resize(mfccs.size(), empty_values);
    // calculate deltas based on 2 neighbouring frames in either direction
    for (int i = 2, last_frame = frames_ - 2; i < last_frame; i++) {
      const ValueList& before_2 = mfccs[i - 2];
      const ValueList& before_1 = mfccs[i - 1];
      const ValueList& after_1 = mfccs[i + 1];
      const ValueList& after_2 = mfccs[i + 2];
      ValueList& delta_frame = delta_mfccs_[i];
      for (int m = 0; m < number_of_mfccs; m++) {
        Value delta = after_1[m] - before_1[m]
                      + 2.0 * (after_2[m] - before_2[m]);
        delta_frame[m] = delta / 10.0;
      }
    }
  }
  return delta_mfccs_;
}

const std::vector<ValueList>& FeatureExtractor::DoubleDeltaMFCCs() {
  if (double_delta_mfccs_.empty()) {
    const std::vector<ValueList>& delta_mfccs = DeltaMFCCs();
    // resize buffers
    const std::size_t number_of_mfccs = delta_mfccs[0].size();
    ValueList empty_values(number_of_mfccs, 0);
    double_delta_mfccs_.resize(delta_mfccs.size(), empty_values);
    // calculate deltas based on 2 neighbouring frames in either direction
    for (int i = 2, last_frame = frames_ - 2; i < last_frame; i++) {
      const ValueList& before_2 = delta_mfccs[i - 2];
      const ValueList& before_1 = delta_mfccs[i - 1];
      const ValueList& after_1 = delta_mfccs[i + 1];
      const ValueList& after_2 = delta_mfccs[i + 2];
      ValueList& double_delta_frame = double_delta_mfccs_[i];
      for (int m = 0; m < number_of_mfccs; m++) {
        Value delta = after_1[m] - before_1[m]
                      + 2.0 * (after_2[m] - before_2[m]);
        double_delta_frame[m] = delta / 10.0;
      }
    }
  }
  return double_delta_mfccs_;
}
  
const std::vector<ValueList>& FeatureExtractor::SpectralFlux() {
  if (spectral_flux_.empty()) {
    // get magnitude spectrum
    const std::vector<ValueList>& spectrum = MagnitudeSpectrum();
    // get magnitude spectrum for each frame
    std::size_t frames = spectrum.size();
    std::size_t bins = spectrum[0].size();
    spectral_flux_.resize(frames);
    // first frame always zeroed
    spectral_flux_[0].assign(bins, 0);
    for (int i = 1; i < frames; i++) {
      const ValueList& spectrum_frame = spectrum[i];
      const ValueList& previous_spectrum_frame = spectrum[i - 1];
      ValueList& flux_frame = spectral_flux_[i];
      flux_frame.resize(bins);
      // get the difference in magnitudes for each bin
      for (int bin = 0; bin < bins; bin++) {
        flux_frame[bin] = spectrum_frame[bin] - previous_spectrum_frame[bin];
      }
      //std::cout << i << ": " << stats::Mean(flux_frame) << '\n';
    }
  }
  return spectral_flux_;
}

} // dsp namespace