Chris@51: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@51: /*
Chris@51:     This file is Copyright (c) 2012 Chris Cannam
Chris@51:   
Chris@51:     Permission is hereby granted, free of charge, to any person
Chris@51:     obtaining a copy of this software and associated documentation
Chris@51:     files (the "Software"), to deal in the Software without
Chris@51:     restriction, including without limitation the rights to use, copy,
Chris@51:     modify, merge, publish, distribute, sublicense, and/or sell copies
Chris@51:     of the Software, and to permit persons to whom the Software is
Chris@51:     furnished to do so, subject to the following conditions:
Chris@51: 
Chris@51:     The above copyright notice and this permission notice shall be
Chris@51:     included in all copies or substantial portions of the Software.
Chris@51: 
Chris@51:     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
Chris@51:     EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
Chris@51:     MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
Chris@51:     NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR
Chris@51:     ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
Chris@51:     CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
Chris@51:     WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
Chris@51: */
Chris@51: 
Chris@73: #ifndef CEPSTRUM_H
Chris@73: #define CEPSTRUM_H
Chris@51: 
Chris@51: #include "vamp-sdk/FFT.h"
Chris@51: #include <cmath>
Chris@51: 
Chris@51: #include <iostream>
Chris@51: #include <exception>
Chris@51: 
Chris@51: class Cepstrum
Chris@51: {
Chris@51: public:
Chris@51:     /**
Chris@51:      * Construct a cepstrum converter based on an n-point FFT.
Chris@51:      */
Chris@51:     Cepstrum(int n) : m_n(n) {
Chris@51: 	if (n & (n-1)) {
Chris@51: 	    throw "N must be a power of two";
Chris@51: 	}
Chris@51:     }
Chris@51:     ~Cepstrum() { }
Chris@51:     
Chris@51:     /**
Chris@51:      * Convert the given frequency-domain data to the cepstral domain.
Chris@51:      *
Chris@51:      * The input must be in the format used for FrequencyDomain data
Chris@51:      * in the Vamp SDK: n/2+1 consecutive pairs of real and imaginary
Chris@51:      * component floats corresponding to bins 0..(n/2) of the FFT
Chris@51:      * output. Thus, n+2 values in total.
Chris@51:      *
Chris@51:      * The output consists of the raw cepstrum of length n.
Chris@51:      *
Chris@51:      * The cepstrum is calculated as the inverse FFT of a
Chris@51:      * synthetically symmetrical base-10 log magnitude spectrum.
Chris@51:      *
Chris@51:      * Returns the mean magnitude of the input spectrum.
Chris@51:      */
Chris@51:     double process(const float *in, double *out) {
Chris@51: 
Chris@51: 	int hs = m_n/2 + 1;
Chris@51: 	double *io = new double[m_n];
Chris@51: 	double *logmag = new double[m_n];
Chris@51: 	double epsilon = 1e-10;
Chris@51: 
Chris@51: 	double magmean = 0.0;
Chris@51: 
Chris@51: 	for (int i = 0; i < hs; ++i) {
Chris@51: 
Chris@73:             double re = in[i*2];
Chris@73:             double im = in[i*2+1];
Chris@73:             double power = re * re + im * im;
Chris@51: 	    double mag = sqrt(power);
Chris@51: 	    magmean += mag;
Chris@51: 
Chris@51: 	    logmag[i] = log10(mag + epsilon);
Chris@51: 
Chris@51: 	    if (i > 0) {
Chris@51: 		// make the log magnitude spectrum symmetrical
Chris@51: 		logmag[m_n - i] = logmag[i];
Chris@51: 	    }
Chris@51: 	}
Chris@51: 	
Chris@51: 	magmean /= hs;
Chris@51: 	/*
Chris@51: 	std::cerr << "logmags:" << std::endl;
Chris@51: 	for (int i = 0; i < m_n; ++i) {
Chris@51: 	    std::cerr << logmag[i] << " ";
Chris@51: 	}
Chris@51: 	std::cerr << std::endl;
Chris@51: 	*/
Chris@51: 	Vamp::FFT::inverse(m_n, logmag, 0, out, io);
Chris@51:     
Chris@51: 	delete[] logmag;
Chris@51: 	delete[] io;
Chris@51: 
Chris@51: 	return magmean;
Chris@51:     }
Chris@51: 
Chris@51: private:
Chris@51:     int m_n;
Chris@51: };
Chris@51: 
Chris@51: #endif