Chris@0: /*
Chris@0:   copyright (C) 2011 I. Irigaray, M. Rocamora
Chris@0: 
Chris@0:   This program is free software: you can redistribute it and/or modify
Chris@0:   it under the terms of the GNU General Public License as published by
Chris@0:   the Free Software Foundation, either version 3 of the License, or
Chris@0:   (at your option) any later version.
Chris@0: 
Chris@0:   This program is distributed in the hope that it will be useful,
Chris@0:   but WITHOUT ANY WARRANTY; without even the implied warranty of
Chris@0:   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Chris@0:   GNU General Public License for more details.
Chris@0: 
Chris@0:   You should have received a copy of the GNU General Public License
Chris@0:   along with this program.  If not, see <http://www.gnu.org/licenses/>.
Chris@0:  */
Chris@0: 
Chris@0: #include "FChTransformF0gram.h"
Chris@0: #include "FChTransformUtils.h"
Chris@0: #include <math.h>
Chris@0: #include <float.h>
Chris@0: //#define DEBUG
Chris@0: #define MAX(x, y) (((x) > (y)) ? (x) : (y))
Chris@0: 
Chris@0: FChTransformF0gram::FChTransformF0gram(float inputSampleRate) :
Chris@0: Plugin(inputSampleRate),
Chris@0: m_currentProgram("default"),
Chris@0: m_stepSize(0), // We are using 0 for step and block size to indicate "not yet set".
Chris@0: m_blockSize(0) {
Chris@0: 
Chris@0:     printf("FUNCTION CALL: %s\n", __FUNCTION__);
Chris@0: 
Chris@0:     m_fs = inputSampleRate;
Chris@0:     // max frequency of interest (Hz)
Chris@0:     m_fmax = 10000.f;
Chris@0:     // warping parameters
Chris@0:     m_warp_params.nsamps_twarp   = 2048;
Chris@0:     //m_warp_params.nsamps_twarp = 8;
Chris@0:     m_warp_params.alpha_max = 4;
Chris@0:     m_warp_params.num_warps = 21;
Chris@0:     //m_warp_params.num_warps = 11;
Chris@0:     m_warp_params.fact_over_samp = 2;
Chris@0:     m_warp_params.alpha_dist = 0;
Chris@0:     // f0 parameters
Chris@0:     m_f0_params.f0min = 80.0;
Chris@0:     m_f0_params.num_octs = 4;
Chris@0:     m_f0_params.num_f0s_per_oct = 192;
Chris@0:     m_f0_params.num_f0_hyps = 5;
Chris@0:     m_f0_params.prefer = true;
Chris@0:     m_f0_params.prefer_mean = 60;
Chris@0:     m_f0_params.prefer_stdev = 18;
Chris@0:     // glogs parameters
Chris@0:     m_glogs_params.HP_logS = true;
Chris@0:     m_glogs_params.att_subharms = 1;
Chris@0: 	// display parameters
Chris@0: 	m_f0gram_mode = true;
Chris@0: 
Chris@0:     m_glogs_params.median_poly_coefs[0] = -0.000000058551680;
Chris@0:     m_glogs_params.median_poly_coefs[1] = -0.000006945207775;
Chris@0:     m_glogs_params.median_poly_coefs[2] = 0.002357223226588;
Chris@0: 
Chris@0:     m_glogs_params.sigma_poly_coefs[0] = 0.000000092782308;
Chris@0:     m_glogs_params.sigma_poly_coefs[1] = 0.000057283574898;
Chris@0:     m_glogs_params.sigma_poly_coefs[2] = 0.022199903714288;
Chris@0: 
Chris@0:     // number of fft points (controls zero-padding)
Chris@0:     m_nfft = m_warp_params.nsamps_twarp;
Chris@0:     // hop in samples
Chris@0:     m_hop = m_warp_params.fact_over_samp * 256;
Chris@0: 
Chris@0:     m_num_f0s = 0;
Chris@0: 
Chris@0: }
Chris@0: 
Chris@0: FChTransformF0gram::~FChTransformF0gram() {
Chris@0:     // remeber to delete everything that deserves to
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: FChTransformF0gram::getIdentifier() const {
Chris@0:     return "fchtransformf0gram";
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: FChTransformF0gram::getName() const {
Chris@0:     return "Fan Chirp Transform F0gram";
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: FChTransformF0gram::getDescription() const {
Chris@0:     // Return something helpful here!
Chris@0:     return "This plug-in produces a representation, called F0gram, which exhibits the salience of the fundamental frequency of the sound sources in the audio file. The computation of the F0gram makes use of the Fan Chirp Transform analysis. It is based on the article \"Fan chirp transform for music representation\"  P. Cancela, E. Lopez, M. Rocamora, International Conference on Digital Audio Effects, 13th. DAFx-10. Graz, Austria - 6-10 Sep 2010.";
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: FChTransformF0gram::getMaker() const {
Chris@0:     // Your name here
Chris@0:     return "Audio Processing Group \n Universidad de la Republica";
Chris@0: }
Chris@0: 
Chris@0: int
Chris@0: FChTransformF0gram::getPluginVersion() const {
Chris@0:     // Increment this each time you release a version that behaves
Chris@0:     // differently from the previous one
Chris@0:     //
Chris@0:     // 0 - initial version from scratch
Chris@0:     return 0;
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: FChTransformF0gram::getCopyright() const {
Chris@0:     // This function is not ideally named.  It does not necessarily
Chris@0:     // need to say who made the plugin -- getMaker does that -- but it
Chris@0:     // should indicate the terms under which it is distributed.  For
Chris@0:     // example, "Copyright (year). All Rights Reserved", or "GPL"
Chris@0:     return "copyright (C) 2011 GPL - Audio Processing Group, UdelaR";
Chris@0: }
Chris@0: 
Chris@0: FChTransformF0gram::InputDomain
Chris@0: FChTransformF0gram::getInputDomain() const {
Chris@0:     return TimeDomain;
Chris@0: }
Chris@0: 
Chris@0: size_t FChTransformF0gram::getPreferredBlockSize() const {
Chris@0:     return 8192; // 0 means "I can handle any block size"
Chris@0: }
Chris@0: 
Chris@0: size_t
Chris@0: FChTransformF0gram::getPreferredStepSize() const {
Chris@0:     return 256; // 0 means "anything sensible"; in practice this
Chris@0:     // means the same as the block size for TimeDomain
Chris@0:     // plugins, or half of it for FrequencyDomain plugins
Chris@0: }
Chris@0: 
Chris@0: size_t
Chris@0: FChTransformF0gram::getMinChannelCount() const {
Chris@0:     return 1;
Chris@0: }
Chris@0: 
Chris@0: size_t
Chris@0: FChTransformF0gram::getMaxChannelCount() const {
Chris@0:     return 1;
Chris@0: }
Chris@0: 
Chris@0: FChTransformF0gram::ParameterList
Chris@0: FChTransformF0gram::getParameterDescriptors() const {
Chris@0:     ParameterList list;
Chris@0: 
Chris@0:     // If the plugin has no adjustable parameters, return an empty
Chris@0:     // list here (and there's no need to provide implementations of
Chris@0:     // getParameter and setParameter in that case either).
Chris@0: 
Chris@0:     // Note that it is your responsibility to make sure the parameters
Chris@0:     // start off having their default values (e.g. in the constructor
Chris@0:     // above).  The host needs to know the default value so it can do
Chris@0:     // things like provide a "reset to default" function, but it will
Chris@0:     // not explicitly set your parameters to their defaults for you if
Chris@0:     // they have not changed in the mean time.
Chris@0: 
Chris@0:     // =============  WARPING PARAMETERS  =============
Chris@0: 
Chris@0:     ParameterDescriptor fmax;
Chris@0:     fmax.identifier = "fmax";
Chris@0:     fmax.name = "Maximum frequency";
Chris@0:     fmax.description = "Maximum frequency of interest for the analysis.";
Chris@0:     fmax.unit = "Hz";
Chris@0:     fmax.minValue = 2000;
Chris@0:     fmax.maxValue = 22050;
Chris@0:     fmax.defaultValue = 10000;
Chris@0:     fmax.isQuantized = true;
Chris@0:     fmax.quantizeStep = 1.0;
Chris@0:     list.push_back(fmax);
Chris@0: 
Chris@0:     ParameterDescriptor nsamp;
Chris@0:     nsamp.identifier = "nsamp";
Chris@0:     nsamp.name = "Number of samples";
Chris@0:     nsamp.description = "Number of samples of the time warped frame";
Chris@0:     nsamp.unit = "samples";
Chris@0:     nsamp.minValue = 128;
Chris@0:     nsamp.maxValue = 4096;
Chris@0:     nsamp.defaultValue = 2048;
Chris@0:     nsamp.isQuantized = true;
Chris@0:     nsamp.quantizeStep = 1.0;
Chris@0:     list.push_back(nsamp);
Chris@0: 
Chris@0:     ParameterDescriptor nfft;
Chris@0:     nfft.identifier = "nfft";
Chris@0:     nfft.name = "FFT number of points";
Chris@0:     nfft.description = "Number of FFT points (controls zero-padding)";
Chris@0:     nfft.unit = "samples";
Chris@0:     nfft.minValue = 0;
Chris@0:     nfft.maxValue = 4;
Chris@0:     nfft.defaultValue = 3;
Chris@0:     nfft.isQuantized = true;
Chris@0:     nfft.quantizeStep = 1.0;
Chris@0:     nfft.valueNames.push_back("256");
Chris@0:     nfft.valueNames.push_back("512");
Chris@0:     nfft.valueNames.push_back("1024");
Chris@0:     nfft.valueNames.push_back("2048");
Chris@0:     nfft.valueNames.push_back("4096");
Chris@0:     nfft.valueNames.push_back("8192");
Chris@0:     list.push_back(nfft);
Chris@0: 
Chris@0:     ParameterDescriptor alpha_max;
Chris@0:     alpha_max.identifier = "alpha_max";
Chris@0:     alpha_max.name = "Maximum alpha value";
Chris@0:     alpha_max.description = "Maximum value for the alpha parameter of the transform.";
Chris@0:     alpha_max.unit = "Hz/s";
Chris@0:     alpha_max.minValue = -10;
Chris@0:     alpha_max.maxValue = 10;
Chris@0:     alpha_max.defaultValue = 5;
Chris@0:     alpha_max.isQuantized = true;
Chris@0:     alpha_max.quantizeStep = 1.0;
Chris@0:     list.push_back(alpha_max);
Chris@0: 
Chris@0:     ParameterDescriptor num_warps;
Chris@0:     num_warps.identifier = "num_warps";
Chris@0:     num_warps.name = "Number of warpings";
Chris@0:     num_warps.description = "Number of different warpings in the specified range (must be odd).";
Chris@0:     num_warps.unit = "";
Chris@0:     num_warps.minValue = 1;
Chris@0:     num_warps.maxValue = 101;
Chris@0:     num_warps.defaultValue = 21;
Chris@0:     num_warps.isQuantized = true;
Chris@0:     num_warps.quantizeStep = 2.0;
Chris@0:     list.push_back(num_warps);
Chris@0: 
Chris@0:     ParameterDescriptor alpha_dist;
Chris@0:     alpha_dist.identifier = "alpha_dist";
Chris@0:     alpha_dist.name = "alpha distribution";
Chris@0:     alpha_dist.description = "Type of distribution of alpha values (linear or log).";
Chris@0:     alpha_dist.unit = "";
Chris@0:     alpha_dist.minValue = 0;
Chris@0:     alpha_dist.maxValue = 1;
Chris@0:     alpha_dist.defaultValue = 1;
Chris@0:     alpha_dist.isQuantized = true;
Chris@0:     alpha_dist.quantizeStep = 1.0;
Chris@0:     // lin (0), log (1)
Chris@0:     alpha_dist.valueNames.push_back("lin");
Chris@0:     alpha_dist.valueNames.push_back("log");
Chris@0:     list.push_back(alpha_dist);
Chris@0: 
Chris@0:     // =============   F0-GRAM PARAMETERS  =============
Chris@0: 
Chris@0:     ParameterDescriptor f0min;
Chris@0:     f0min.identifier = "f0min";
Chris@0:     f0min.name = "min f0";
Chris@0:     f0min.description = "Minimum fundamental frequency (f0) value.";
Chris@0:     f0min.unit = "Hz";
Chris@0:     f0min.minValue = 1;
Chris@0:     f0min.maxValue = 500;
Chris@0:     f0min.defaultValue = 80;
Chris@0:     f0min.isQuantized = true;
Chris@0:     f0min.quantizeStep = 1.0;
Chris@0:     list.push_back(f0min);
Chris@0: 
Chris@0:     ParameterDescriptor num_octs;
Chris@0:     num_octs.identifier = "num_octs";
Chris@0:     num_octs.name = "number of octaves";
Chris@0:     num_octs.description = "Number of octaves for F0gram computation.";
Chris@0:     num_octs.unit = "";
Chris@0:     num_octs.minValue = 1;
Chris@0:     num_octs.maxValue = 10;
Chris@0:     num_octs.defaultValue = 4;
Chris@0:     num_octs.isQuantized = true;
Chris@0:     num_octs.quantizeStep = 1.0;
Chris@0:     list.push_back(num_octs);
Chris@0: 
Chris@0:     ParameterDescriptor num_f0_hyps;
Chris@0:     num_f0_hyps.identifier = "num_f0_hyps";
Chris@0:     num_f0_hyps.name = "number of f0 hypotesis";
Chris@0:     num_f0_hyps.description = "Number of f0 hypotesis to extract.";
Chris@0:     num_f0_hyps.unit = "";
Chris@0:     num_f0_hyps.minValue = 1;
Chris@0:     num_f0_hyps.maxValue = 100;
Chris@0:     num_f0_hyps.defaultValue = 10;
Chris@0:     num_f0_hyps.isQuantized = true;
Chris@0:     num_f0_hyps.quantizeStep = 1.0;
Chris@0:     list.push_back(num_f0_hyps);
Chris@0: 
Chris@0:     ParameterDescriptor f0s_per_oct;
Chris@0:     f0s_per_oct.identifier = "f0s_per_oct";
Chris@0:     f0s_per_oct.name = "f0 values per octave";
Chris@0:     f0s_per_oct.description = "Number of f0 values per octave.";
Chris@0:     f0s_per_oct.unit = "";
Chris@0:     f0s_per_oct.minValue = 12;
Chris@0:     f0s_per_oct.maxValue = 768;
Chris@0:     f0s_per_oct.defaultValue = 192;
Chris@0:     f0s_per_oct.isQuantized = true;
Chris@0:     f0s_per_oct.quantizeStep = 1.0;
Chris@0:     list.push_back(f0s_per_oct);
Chris@0: 
Chris@0:     ParameterDescriptor f0_prefer_fun;
Chris@0:     f0_prefer_fun.identifier = "f0_prefer_fun";
Chris@0:     f0_prefer_fun.name = "f0 preference function";
Chris@0:     f0_prefer_fun.description = "Whether to use a f0 weighting function.";
Chris@0:     f0_prefer_fun.unit = "";
Chris@0:     f0_prefer_fun.minValue = 0;
Chris@0:     f0_prefer_fun.maxValue = 1;
Chris@0:     f0_prefer_fun.defaultValue = 1;
Chris@0:     f0_prefer_fun.isQuantized = true;
Chris@0:     f0_prefer_fun.quantizeStep = 1.0;
Chris@0:     list.push_back(f0_prefer_fun);
Chris@0: 
Chris@0:     ParameterDescriptor f0_prefer_mean;
Chris@0:     f0_prefer_mean.identifier = "f0_prefer_mean";
Chris@0:     f0_prefer_mean.name = "mean f0 preference function";
Chris@0:     f0_prefer_mean.description = "Mean value for f0 weighting function (MIDI number).";
Chris@0:     f0_prefer_mean.unit = "";
Chris@0:     f0_prefer_mean.minValue = 1;
Chris@0:     f0_prefer_mean.maxValue = 127;
Chris@0:     f0_prefer_mean.defaultValue = 60;
Chris@0:     f0_prefer_mean.isQuantized = true;
Chris@0:     f0_prefer_mean.quantizeStep = 1.0;
Chris@0:     list.push_back(f0_prefer_mean);
Chris@0: 
Chris@0:     ParameterDescriptor f0_prefer_stdev;
Chris@0:     f0_prefer_stdev.identifier = "f0_prefer_stdev";
Chris@0:     f0_prefer_stdev.name = "stdev of f0 preference function";
Chris@0:     f0_prefer_stdev.description = "Stdev for f0 weighting function (MIDI number).";
Chris@0:     f0_prefer_stdev.unit = "";
Chris@0:     f0_prefer_stdev.minValue = 1;
Chris@0:     f0_prefer_stdev.maxValue = 127;
Chris@0:     f0_prefer_stdev.defaultValue = 18;
Chris@0:     f0_prefer_stdev.isQuantized = true;
Chris@0:     f0_prefer_stdev.quantizeStep = 1.0;
Chris@0:     list.push_back(f0_prefer_stdev);
Chris@0: 
Chris@0:     ParameterDescriptor f0gram_mode;
Chris@0:     f0gram_mode.identifier = "f0gram_mode";
Chris@0:     f0gram_mode.name = "display mode of f0gram";
Chris@0:     f0gram_mode.description = "Display all bins of the best direction, or the best bin for each direction.";
Chris@0:     f0gram_mode.unit = "";
Chris@0:     f0gram_mode.minValue = 0;
Chris@0:     f0gram_mode.maxValue = 1;
Chris@0:     f0gram_mode.defaultValue = 1;
Chris@0:     f0gram_mode.isQuantized = true;
Chris@0:     f0gram_mode.quantizeStep = 1.0;
Chris@0:     list.push_back(f0gram_mode);
Chris@0: 
Chris@0:     return list;
Chris@0: }
Chris@0: 
Chris@0: float
Chris@0: FChTransformF0gram::getParameter(string identifier) const {
Chris@0: 
Chris@0:     printf("FUNCTION CALL: %s(%s)\n", __FUNCTION__, identifier.c_str());
Chris@0: 
Chris@0:     if (identifier == "fmax") {
Chris@0:         return m_fmax;
Chris@0:     } else if (identifier == "nsamp") {
Chris@0:         return m_warp_params.nsamps_twarp;
Chris@0:     } else if (identifier == "alpha_max") {
Chris@0:         return m_warp_params.alpha_max;
Chris@0:     } else if (identifier == "num_warps") {
Chris@0:         return m_warp_params.num_warps;
Chris@0:     } else if (identifier == "alpha_dist") {
Chris@0:         return m_warp_params.alpha_dist;
Chris@0:     } else if (identifier == "nfft") {
Chris@0:         return m_nfft;
Chris@0:     } else if (identifier == "f0min") {
Chris@0:         return m_f0_params.f0min;
Chris@0:     } else if (identifier == "num_octs") {
Chris@0:         return m_f0_params.num_octs;
Chris@0:     } else if (identifier == "f0s_per_oct") {
Chris@0:         return m_f0_params.num_f0s_per_oct;
Chris@0:     } else if (identifier == "num_f0_hyps") {
Chris@0:         return m_f0_params.num_f0_hyps;
Chris@0:     } else if (identifier == "f0_prefer_fun") {
Chris@0:         return m_f0_params.prefer;
Chris@0:     } else if (identifier == "f0_prefer_mean") {
Chris@0:         return m_f0_params.prefer_mean;
Chris@0:     } else if (identifier == "f0_prefer_stdev") {
Chris@0:         return m_f0_params.prefer_stdev;
Chris@0: 	} else if (identifier == "f0gram_mode") {
Chris@0:         return m_f0gram_mode;
Chris@0:     } else {
Chris@0:         return 0.f;
Chris@0:     }
Chris@0: 
Chris@0: }
Chris@0: 
Chris@0: void FChTransformF0gram::setParameter(string identifier, float value) {
Chris@0: 
Chris@0:     printf("FUNCTION CALL: %s(%s) = %f.\n", __FUNCTION__, identifier.c_str(), value);
Chris@0: 
Chris@0:     if (identifier == "fmax") {
Chris@0:         m_fmax = value;
Chris@0:     } else if (identifier == "nsamp") {
Chris@0:         m_warp_params.nsamps_twarp = value;
Chris@0:     } else if (identifier == "alpha_max") {
Chris@0:         m_warp_params.alpha_max = value;
Chris@0:     } else if (identifier == "num_warps") {
Chris@0:         m_warp_params.num_warps = value;
Chris@0:     } else if (identifier == "alpha_dist") {
Chris@0:         m_warp_params.alpha_dist = value;
Chris@0:     } else if (identifier == "nfft") {
Chris@0:         m_nfft = value;
Chris@0:     } else if (identifier == "f0min") {
Chris@0:         m_f0_params.f0min = value;
Chris@0:     } else if (identifier == "num_octs") {
Chris@0:         m_f0_params.num_octs = value;
Chris@0:     } else if (identifier == "f0s_per_oct") {
Chris@0:         m_f0_params.num_f0s_per_oct = value;
Chris@0:     } else if (identifier == "num_f0_hyps") {
Chris@0:         m_f0_params.num_f0_hyps = value;
Chris@0:     } else if (identifier == "f0_prefer_fun") {
Chris@0:         m_f0_params.prefer = value;
Chris@0:     } else if (identifier == "f0_prefer_mean") {
Chris@0:         m_f0_params.prefer_mean = value;
Chris@0:     } else if (identifier == "f0_prefer_stdev") {
Chris@0:         m_f0_params.prefer_stdev = value;
Chris@0:     } else if (identifier == "f0gram_mode") {
Chris@0:         m_f0gram_mode = value;
Chris@0:     }
Chris@0: 
Chris@0: }
Chris@0: 
Chris@0: FChTransformF0gram::ProgramList
Chris@0: FChTransformF0gram::getPrograms() const {
Chris@0:     ProgramList list;
Chris@0: 
Chris@0:     list.push_back("default");
Chris@0: 
Chris@0:     return list;
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: FChTransformF0gram::getCurrentProgram() const {
Chris@0:     return m_currentProgram;
Chris@0: }
Chris@0: 
Chris@0: void
Chris@0: FChTransformF0gram::selectProgram(string name) {
Chris@0: 
Chris@0:     printf("FUNCTION CALL: %s\n", __FUNCTION__);
Chris@0: 
Chris@0:     m_currentProgram = name;
Chris@0: 
Chris@0:     if (name == "default") {
Chris@0:         m_fmax = 10000.f;
Chris@0: 
Chris@0:         m_warp_params.nsamps_twarp = 2048;
Chris@0:         m_warp_params.alpha_max = 4;
Chris@0:         m_warp_params.num_warps = 21;
Chris@0:         m_warp_params.fact_over_samp = 2;
Chris@0:         m_warp_params.alpha_dist = 0;
Chris@0: 
Chris@0:         m_f0_params.f0min = 80.0;
Chris@0:         m_f0_params.num_octs = 4;
Chris@0:         m_f0_params.num_f0s_per_oct = 192;
Chris@0:         m_f0_params.num_f0_hyps = 5;
Chris@0:         m_f0_params.prefer = true;
Chris@0:         m_f0_params.prefer_mean = 60;
Chris@0:         m_f0_params.prefer_stdev = 18;
Chris@0: 
Chris@0:         m_glogs_params.HP_logS = true;
Chris@0:         m_glogs_params.att_subharms = 1;
Chris@0: 
Chris@0:         m_glogs_params.median_poly_coefs[0] = -0.000000058551680;
Chris@0:         m_glogs_params.median_poly_coefs[1] = -0.000006945207775;
Chris@0:         m_glogs_params.median_poly_coefs[2] = 0.002357223226588;
Chris@0: 
Chris@0:         m_glogs_params.sigma_poly_coefs[0] = 0.000000092782308;
Chris@0:         m_glogs_params.sigma_poly_coefs[1] = 0.000057283574898;
Chris@0:         m_glogs_params.sigma_poly_coefs[2] = 0.022199903714288;
Chris@0: 
Chris@0:         m_nfft = m_warp_params.nsamps_twarp;
Chris@0:         m_hop = m_warp_params.fact_over_samp * 256;
Chris@0: 
Chris@0:         m_num_f0s = 0;
Chris@0: 
Chris@0: 		m_f0gram_mode = 1;
Chris@0: 
Chris@0:     }
Chris@0: }
Chris@0: 
Chris@0: FChTransformF0gram::OutputList
Chris@0: FChTransformF0gram::getOutputDescriptors() const {
Chris@0: 
Chris@0:     printf("FUNCTION CALL: %s\n", __FUNCTION__);
Chris@0: 
Chris@0:     OutputList list;
Chris@0: 
Chris@0:     // See OutputDescriptor documentation for the possibilities here.
Chris@0:     // Every plugin must have at least one output.
Chris@0: 
Chris@0:     /* f0 values of F0gram grid as string values */
Chris@0:     vector<string> f0values;
Chris@0:     size_t ind = 0;
Chris@0:     char f0String[10];
Chris@0:     while (ind < m_num_f0s) {
Chris@0:         sprintf(f0String, "%4.2f", m_f0s[ind]);
Chris@0:         f0values.push_back(f0String);
Chris@0:         ind++;
Chris@0:     }
Chris@0: 
Chris@0:     /* The F0gram */
Chris@0:     OutputDescriptor d;
Chris@0:     d.identifier = "f0gram";
Chris@0:     d.name = "F0gram: salience of f0s";
Chris@0:     d.description = "This representation show the salience of the different f0s in the signal.";
Chris@0:     d.unit = "Hertz";
Chris@0:     d.hasFixedBinCount = true;
Chris@0:     //d.binCount = m_num_f0s;
Chris@0: 	//d.binCount = m_blockSize/2+1;
Chris@0: 	//d.binCount = m_warp_params.nsamps_twarp/2+1;
Chris@0: 	//d.binCount = m_warpings.nsamps_torig;
Chris@0: 	d.binCount = m_f0_params.num_octs*m_f0_params.num_f0s_per_oct;
Chris@0:     d.binNames = f0values;
Chris@0:     d.hasKnownExtents = false;
Chris@0:     d.isQuantized = false;
Chris@0:     d.sampleType = OutputDescriptor::OneSamplePerStep;
Chris@0:     d.hasDuration = false;
Chris@0:     list.push_back(d);
Chris@0: 
Chris@0:     return list;
Chris@0: }
Chris@0: 
Chris@0: bool
Chris@0: FChTransformF0gram::initialise(size_t channels, size_t stepSize, size_t blockSize) {
Chris@0:     if (channels < getMinChannelCount() ||
Chris@0:             channels > getMaxChannelCount()) return false;
Chris@0: 
Chris@0:     printf("FUNCTION CALL: %s\n", __FUNCTION__);
Chris@0: 
Chris@0:     // set blockSize and stepSize (but changed below)
Chris@0:     m_blockSize = blockSize;
Chris@0:     m_stepSize = stepSize;
Chris@0: 
Chris@0:     // WARNING !!!
Chris@0:     // these values in fact are determined by the sampling frequency m_fs
Chris@0:     // the parameters used below correspond to default values i.e. m_fs = 44.100 Hz
Chris@0:     //m_blockSize = 4 * m_warp_params.nsamps_twarp;
Chris@0:     m_stepSize = floor(m_hop / m_warp_params.fact_over_samp);
Chris@0: 
Chris@0:     /* initialise m_warp_params  */
Chris@0:     //    FChTF0gram:warping_design m_warpings = new warping_design;
Chris@0:     /* initialise m_f0_params    */
Chris@0: 
Chris@0:     /* initialise m_glogs_params */
Chris@0: 	design_GLogS();
Chris@0: 
Chris@0:     /* design of FChT */
Chris@0:     // design_fcht(m_warps, m_accums, m_f0s)
Chris@0:     design_FChT();
Chris@0: 
Chris@0: 	design_FFT();
Chris@0: 
Chris@0: 	design_LPF();
Chris@0: 
Chris@0: 	design_time_window();
Chris@0: 
Chris@0: 	// Create Hanning window for warped signals
Chris@0: 	mp_HanningWindow = new double[m_warp_params.nsamps_twarp];
Chris@0: 	bool normalize = false;
Chris@0: 	hanning_window(mp_HanningWindow, m_warp_params.nsamps_twarp, normalize);
Chris@0: 
Chris@0: 	printf("	End of initialise().\n");
Chris@0:     return true;
Chris@0: }
Chris@0: 
Chris@0: void
Chris@0: FChTransformF0gram::design_GLogS() {
Chris@0: 	printf("	Running design_GLogS().\n");
Chris@0: 
Chris@0: 	// total number & initial quantity of f0s
Chris@0: 	m_glogs_init_f0s = (size_t)(((double)m_f0_params.num_f0s_per_oct)*log2(5.0))+1;
Chris@0: 	m_glogs_num_f0s = (m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct + m_glogs_init_f0s;
Chris@0: 
Chris@0: 	// Initialize arrays
Chris@0: 	m_glogs_f0 = new double[m_glogs_num_f0s];
Chris@0: 	m_glogs = new double[m_glogs_num_f0s*m_warp_params.num_warps];
Chris@0: 	m_glogs_n = new size_t[m_glogs_num_f0s];
Chris@0: 	m_glogs_index = new size_t[m_glogs_num_f0s];
Chris@0: 
Chris@0: 	// Compute f0 values
Chris@0: 	m_glogs_harmonic_count = 0;
Chris@0: 	double factor = (double)(m_warp_params.nsamps_twarp/2)/(double)(m_warp_params.nsamps_twarp/2+1);
Chris@0: 	for (size_t i = 0; i < m_glogs_num_f0s; i++) {
Chris@0: 		m_glogs_f0[i] = (m_f0_params.f0min/5.0)*pow(2.0,(double)i/(double)m_f0_params.num_f0s_per_oct);
Chris@0: 		// for every f0 compute number of partials less or equal than m_fmax.
Chris@0: 		m_glogs_n[i] = m_fmax*factor/m_glogs_f0[i];
Chris@0: 		m_glogs_index[i] = m_glogs_harmonic_count; 
Chris@0: 		m_glogs_harmonic_count += m_glogs_n[i];
Chris@0: 	}
Chris@0: 
Chris@0: 	// Initialize arrays for interpolation
Chris@0: 	m_glogs_posint = new size_t[m_glogs_harmonic_count];
Chris@0: 	m_glogs_posfrac = new double[m_glogs_harmonic_count];
Chris@0: 	m_glogs_interp = new double[m_glogs_harmonic_count];
Chris@0: 
Chris@0: 	// Compute int & frac of interpolation positions
Chris@0: 	size_t aux_index = 0;
Chris@0: 	double aux_pos;
Chris@0: 	for (size_t i = 0; i < m_glogs_num_f0s; i++) {
Chris@0: 		for (size_t j = 1; j <= m_glogs_n[i]; j++) {
Chris@0: 			// indice en el vector de largo t_warp/2+1 donde el ultimo valor corresponde a f=m_fmax
Chris@0: 			aux_pos = ((double)j*m_glogs_f0[i])*((double)(m_warp_params.nsamps_twarp/2+1))/m_fmax;
Chris@0: 			m_glogs_posint[aux_index] = (size_t)aux_pos;
Chris@0: 			m_glogs_posfrac[aux_index] = aux_pos - (double)m_glogs_posint[aux_index];
Chris@0: 			aux_index++;
Chris@0: 		}
Chris@0: 	}
Chris@0: 
Chris@0: 	// Third harmonic attenuation
Chris@0: 	double aux_third_harmonic;
Chris@0: 	m_glogs_third_harmonic_posint = new size_t[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct];
Chris@0: 	m_glogs_third_harmonic_posfrac = new double[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct];
Chris@0: 	for (size_t i = 0; i < (m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct; i++) {
Chris@0: 		aux_third_harmonic = (double)i + (double)m_glogs_init_f0s - ((double)m_f0_params.num_f0s_per_oct)*log2(3.0);
Chris@0: 			//printf("	aux_third_harmonic = %f.\n", aux_third_harmonic);
Chris@0: 		m_glogs_third_harmonic_posint[i] = (size_t)aux_third_harmonic;
Chris@0: 		m_glogs_third_harmonic_posfrac[i] = aux_third_harmonic - (double)(m_glogs_third_harmonic_posint[i]);
Chris@0: 	}
Chris@0: 	m_glogs_third_harmonic = new double[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct];
Chris@0: 
Chris@0: 	// Fifth harmonic attenuation
Chris@0: 	double aux_fifth_harmonic;
Chris@0: 	m_glogs_fifth_harmonic_posint = new size_t[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct];
Chris@0: 	m_glogs_fifth_harmonic_posfrac = new double[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct];
Chris@0: 	for (size_t i = 0; i < (m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct; i++) {
Chris@0: 		aux_fifth_harmonic = (double)i + (double)m_glogs_init_f0s - ((double)m_f0_params.num_f0s_per_oct)*log2(5.0);
Chris@0: 			//printf("	aux_fifth_harmonic = %f.\n", aux_fifth_harmonic);
Chris@0: 		m_glogs_fifth_harmonic_posint[i] = (size_t)aux_fifth_harmonic;
Chris@0: 		m_glogs_fifth_harmonic_posfrac[i] = aux_fifth_harmonic - (double)(m_glogs_fifth_harmonic_posint[i]);
Chris@0: 	}
Chris@0: 	m_glogs_fifth_harmonic = new double[(m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct];
Chris@0: 
Chris@0: 	// Normalization & attenuation windows
Chris@0: 	m_glogs_f0_preference_weights = new double[m_f0_params.num_octs*m_f0_params.num_f0s_per_oct];
Chris@0: 	m_glogs_median_correction = new double[m_f0_params.num_octs*m_f0_params.num_f0s_per_oct];
Chris@0: 	m_glogs_sigma_correction = new double[m_f0_params.num_octs*m_f0_params.num_f0s_per_oct];
Chris@0: 	m_glogs_hf_smoothing_window = new double[m_warp_params.nsamps_twarp/2+1];
Chris@0: 	double MIDI_value;
Chris@0: 	for (size_t i = 0; i < m_f0_params.num_octs*m_f0_params.num_f0s_per_oct; i++) {
Chris@0: 		MIDI_value = 69.0 + 12.0 * log2(m_glogs_f0[i + m_glogs_init_f0s]/440.0);
Chris@0: 		m_glogs_f0_preference_weights[i] = 1.0/sqrt(2.0*M_PI*m_f0_params.prefer_stdev*m_f0_params.prefer_stdev)*exp(-(MIDI_value-m_f0_params.prefer_mean)*(MIDI_value-m_f0_params.prefer_mean)/(2.0*m_f0_params.prefer_stdev*m_f0_params.prefer_stdev));
Chris@0: 		m_glogs_f0_preference_weights[i] = (0.01 + m_glogs_f0_preference_weights[i]) / (1.01);
Chris@0: 		
Chris@0: 		m_glogs_median_correction[i] = m_glogs_params.median_poly_coefs[0]*(i+1.0)*(i+1.0) + m_glogs_params.median_poly_coefs[1]*(i+1.0) + m_glogs_params.median_poly_coefs[2];
Chris@0: 		m_glogs_sigma_correction[i] = 1.0 / (m_glogs_params.sigma_poly_coefs[0]*(i+1.0)*(i+1.0) + m_glogs_params.sigma_poly_coefs[1]*(i+1.0) + m_glogs_params.sigma_poly_coefs[2]);
Chris@0: 	}
Chris@0: 	
Chris@0: 	double smooth_width = 1000.0; // hertz.
Chris@0: 	double smooth_aux = (double)(m_warp_params.nsamps_twarp/2+1)*(m_fmax-smooth_width)/m_fmax;
Chris@0: 	for (size_t i = 0; i < m_warp_params.nsamps_twarp/2+1; i++) {
Chris@0: 		if (i <  smooth_aux) {
Chris@0: 			m_glogs_hf_smoothing_window[i] = 1.0;
Chris@0: 		} else {
Chris@0: 			m_glogs_hf_smoothing_window[i] = ((double)i - (double)m_warp_params.nsamps_twarp/2.0)*(-1.0/((double)(m_warp_params.nsamps_twarp/2+1)-smooth_aux));
Chris@0: 		}
Chris@0: 	}
Chris@0: 
Chris@0: 	printf("	End of design_GLogS().\n");
Chris@0: 	
Chris@0: }
Chris@0: 
Chris@0: void
Chris@0: FChTransformF0gram::design_FFT() {
Chris@0: 	printf("	Running design_FFT().\n");
Chris@0: 
Chris@0:     in = (fftw_complex*) fftw_malloc(sizeof (fftw_complex) * m_nfft);
Chris@0:     out = (fftw_complex*) fftw_malloc(sizeof (fftw_complex) * m_nfft);
Chris@0: 	//TODO verificar que el tipo de datos de in_window es del tipo double, era float.
Chris@0:     in_window = (double*) fftw_malloc(sizeof (double) * m_nfft);
Chris@0:     planFFT = fftw_plan_dft_1d(m_nfft, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
Chris@0: 
Chris@0: 	//TODO hacer diseño del FFT para el filtrado pasabajos.
Chris@0: 
Chris@0: }
Chris@0: 
Chris@0: void
Chris@0: FChTransformF0gram::design_FChT() {
Chris@0: 
Chris@0: 	printf("	Running design_FChT().\n");
Chris@0:     /*
Chris@0:      * FILES FOR DEBUGGING
Chris@0:      */
Chris@0: 
Chris@0:     //ofstream output("output.txt");
Chris@0: 
Chris@0: 
Chris@0:     /*  =============  WARPING DESIGN   ============= */
Chris@0: 
Chris@0:     // sampling frequency after oversampling
Chris@0:     m_warpings.fs_orig = m_warp_params.fact_over_samp * m_fs;
Chris@0: 
Chris@0:     // number of samples of the original signal frame
Chris@0:     m_warpings.nsamps_torig = 4 * m_warp_params.fact_over_samp * m_warp_params.nsamps_twarp;
Chris@0:     // equivalent to: m_warpings.nsamps_torig = m_warp_params.fact_over_samp * m_blockSize;
Chris@0: 
Chris@0:     // time instants of the original signal frame
Chris@0:     double t_orig[m_warpings.nsamps_torig];
Chris@0:     //float * t_orig = new float [m_warpings.nsamps_torig];
Chris@0:     for (size_t ind = 0; ind < m_warpings.nsamps_torig; ind++) {
Chris@0:         t_orig[ind] = ((double)(ind + 1) - (double)m_warpings.nsamps_torig / 2.0) / m_warpings.fs_orig;
Chris@0:     }
Chris@0: 
Chris@0:     // linear chirps warping definition as relative frequency deviation
Chris@0:     //double * freq_relative = new double [m_warpings.nsamps_torig * m_warp_params.num_warps];
Chris@0: 	//TODO
Chris@0: 	double *freq_relative = new double [m_warpings.nsamps_torig * m_warp_params.num_warps];
Chris@0:     define_warps_linear_chirps(freq_relative, t_orig);
Chris@0: 
Chris@0:     // maximum relative frequency deviation
Chris@0:     double freq_relative_max = 0;
Chris@0:     for (size_t i = 0; i < m_warpings.nsamps_torig; i++)
Chris@0:         for (size_t j = 0; j < m_warp_params.num_warps; j++)
Chris@0:             if (freq_relative_max < freq_relative[j * m_warpings.nsamps_torig + i])
Chris@0:                 freq_relative_max = freq_relative[j * m_warpings.nsamps_torig + i];
Chris@0: 
Chris@0:     // sampling frequency of warped signal to be free of aliasing up to fmax
Chris@0:     m_warpings.fs_warp = 2 * m_fmax * freq_relative_max;
Chris@0: 
Chris@0:     // time instants of the warped signal frame
Chris@0:     double t_warp[m_warp_params.nsamps_twarp];
Chris@0:     for (size_t ind = 0; ind < m_warp_params.nsamps_twarp; ind++) {
Chris@0:         t_warp[ind] = ((double)((int)(ind + 1)- (int)m_warp_params.nsamps_twarp / 2)) / (double)m_warpings.fs_warp;
Chris@0:     }
Chris@0: 
Chris@0:     // design of warpings for efficient interpolation
Chris@0:     design_warps(freq_relative, t_orig, t_warp);
Chris@0: 
Chris@0: 
Chris@0:     /*
Chris@0:      * FILES FOR DEBUGGING
Chris@0:      */
Chris@0: 
Chris@0:     /*
Chris@0:     output << "chirp_rates" << endl;
Chris@0:     for (size_t j = 0; j < m_warp_params.num_warps; j++){
Chris@0:         output << m_warpings.chirp_rates[j];
Chris@0:         output << " ";
Chris@0:     }
Chris@0:     output << endl << "freq_relative" << endl;
Chris@0: 
Chris@0:     for (size_t i = 0; i < m_warpings.nsamps_torig; i++){
Chris@0:             for (size_t j = 0; j < m_warp_params.num_warps; j++){
Chris@0:                    output << freq_relative[j * m_warpings.nsamps_torig + i];
Chris@0:                    output << " ";
Chris@0:             }
Chris@0:             output << endl;
Chris@0:     }
Chris@0: 
Chris@0:     output << endl << "t_orig" << endl;
Chris@0: 
Chris@0:     for (size_t i = 0; i < m_warpings.nsamps_torig; i++){
Chris@0:                   output << t_orig[i] << endl ;
Chris@0:     }
Chris@0:      */
Chris@0: 
Chris@0: 	delete [] freq_relative;
Chris@0:     //output.close();
Chris@0: 
Chris@0:     /*  =============  FFTW PLAN DESIGN   ============= */
Chris@0: 	// Initialize 2-d array for warped signals
Chris@0: 	x_warping = new double[m_warp_params.nsamps_twarp];
Chris@0: 	m_absFanChirpTransform = (double*)fftw_malloc(sizeof (double) * m_warp_params.num_warps * (m_warp_params.nsamps_twarp/2 + 1));
Chris@0: 	m_auxFanChirpTransform = (fftw_complex*)fftw_malloc(sizeof ( fftw_complex) * (m_warp_params.nsamps_twarp/2 + 1));
Chris@0: 	plan_forward_xwarping = fftw_plan_dft_r2c_1d(m_warp_params.nsamps_twarp, x_warping, m_auxFanChirpTransform, FFTW_ESTIMATE);
Chris@0: 
Chris@0: 	printf("	End of design_FChT().\n");
Chris@0: 
Chris@0: }
Chris@0: 
Chris@0: void
Chris@0: FChTransformF0gram::design_warps(double * freq_relative, double * t_orig, double * t_warp) {
Chris@0:     /* the warping is done by interpolating the original signal in time instants
Chris@0:        given by the desired frequency deviation, to do this, the interpolation
Chris@0:        instants are stored in a structure as an integer index and a fractional value
Chris@0:        hypothesis: sampling frequency at the central point equals the original
Chris@0:      */
Chris@0: 	printf("	Running design_warps().\n");
Chris@0: 
Chris@0:     m_warpings.pos_int = new size_t[m_warp_params.num_warps * m_warp_params.nsamps_twarp];
Chris@0: 	m_warpings.pos_frac = new double[m_warp_params.num_warps * m_warp_params.nsamps_twarp];
Chris@0: 
Chris@0: 	// vector of phase values
Chris@0: 	double *phi = new double[m_warpings.nsamps_torig];
Chris@0: 	double aux;
Chris@0: 
Chris@0: 	// warped positions
Chris@0: 	double *pos1 = new double[m_warp_params.nsamps_twarp*m_warp_params.num_warps];
Chris@0: 	
Chris@0:     for (size_t i = 0; i < m_warp_params.num_warps; i++) {
Chris@0:         // vector of phase values
Chris@0:         //    float * phi;
Chris@0:         // integration of relative frequency to obtain phase values
Chris@0:         // phi = cumtrapz(t_orig,freq_relative(:,i)');
Chris@0:         // centering of phase values to force original frequency in the middle
Chris@0:         //phi = phi - phi(end/2);
Chris@0:         // interpolation of phase values to obtain warped positions
Chris@0:         //pos1(i,:) = interp1(phi,t_orig,t_warp)*fs_orig + length(t_orig)/2;
Chris@0: 		
Chris@0: 		// integration of relative frequency to obtain phase values
Chris@0: 		cumtrapz(t_orig, freq_relative + i*(m_warpings.nsamps_torig), m_warpings.nsamps_torig, phi);
Chris@0: 
Chris@0: 		// centering of phase values to force original frequency in the middle
Chris@0: 		aux = phi[m_warpings.nsamps_torig/2];
Chris@0: 		for (size_t j = 0; j < m_warpings.nsamps_torig; j++) {
Chris@0: 			phi[j] -= aux;
Chris@0: 		} //for
Chris@0: 
Chris@0: 		// interpolation of phase values to obtain warped positions
Chris@0: 		interp1(phi, t_orig, m_warpings.nsamps_torig, t_warp, pos1 + i*m_warp_params.nsamps_twarp, m_warp_params.nsamps_twarp);
Chris@0: 		
Chris@0:     }
Chris@0: 
Chris@0:     // % previous sample index
Chris@0:     // pos1_int = uint32(floor(pos1))';
Chris@0:     // % integer corresponding to previous sample index in "c"
Chris@0:     // warps.pos1_int = (pos1_int - uint32(1));
Chris@0:     // % fractional value that defines the warped position
Chris@0:     // warps.pos1_frac = (double(pos1)' - double(pos1_int));
Chris@0: 
Chris@0: 	// m_warpings.pos_int = new size_t[m_warp_params.num_warps * m_warp_params.nsamps_twarp];
Chris@0: 	for (size_t j = 0; j < m_warp_params.nsamps_twarp*m_warp_params.num_warps; j++) {
Chris@0: 		// previous sample index
Chris@0: 		pos1[j] = pos1[j]*m_warpings.fs_orig + m_warpings.nsamps_torig/2 + 1;
Chris@0: 		m_warpings.pos_int[j] = (size_t) pos1[j];
Chris@0: 		m_warpings.pos_frac[j] = pos1[j] - (double)(m_warpings.pos_int[j]);
Chris@0: 	} //for
Chris@0: 
Chris@0: 	delete [] phi;
Chris@0: 	delete [] pos1;
Chris@0: 
Chris@0: 	printf("	End of design_warps().\n");
Chris@0: 
Chris@0: }
Chris@0: 
Chris@0: void
Chris@0: FChTransformF0gram::define_warps_linear_chirps(double * freq_relative, double * t_orig) {
Chris@0:     /**  define warps as relative frequency deviation from original frequency
Chris@0:                  t_orig : time vector
Chris@0:           freq_relative : relative frequency deviations
Chris@0:      */
Chris@0: 	printf("	Running define_warps_linear_chirps().\n");
Chris@0: 
Chris@0:     if (m_warp_params.alpha_dist == 0) {
Chris@0: 
Chris@0:         // linear alpha values spacing
Chris@0:         m_warpings.chirp_rates = new double [m_warp_params.num_warps];
Chris@0:         // WARNING m_warp_params.num_warps must be odd
Chris@0:         m_warpings.chirp_rates[0] = -m_warp_params.alpha_max;
Chris@0:         double increment = (double) m_warp_params.alpha_max / ((m_warp_params.num_warps - 1) / 2);
Chris@0: 
Chris@0:         for (size_t ind = 1; ind < m_warp_params.num_warps; ind++) {
Chris@0:             m_warpings.chirp_rates[ind] = m_warpings.chirp_rates[ind - 1] + increment;
Chris@0:         }
Chris@0:         // force zero value
Chris@0:         m_warpings.chirp_rates[(int) ((m_warp_params.num_warps - 1) / 2)] = 0;
Chris@0: 
Chris@0:     } else {
Chris@0:         // log alpha values spacing
Chris@0:         m_warpings.chirp_rates = new double [m_warp_params.num_warps];
Chris@0: 
Chris@0:         // force zero value
Chris@0:         int middle_point = (int) ((m_warp_params.num_warps - 1) / 2);
Chris@0:         m_warpings.chirp_rates[middle_point] = 0;
Chris@0: 
Chris@0:         double logMax = log10(m_warp_params.alpha_max + 1);
Chris@0:         double increment = logMax / ((m_warp_params.num_warps - 1) / 2.0f);
Chris@0:         double exponent = 0;
Chris@0: 
Chris@0:         // fill positive values
Chris@0:         int ind_log = middle_point;
Chris@0:         for (size_t ind = 0; ind < (m_warp_params.num_warps + 1) / 2; ind++) {
Chris@0:             m_warpings.chirp_rates[ind_log] = pow(10, exponent) - 1;
Chris@0:             exponent += increment;
Chris@0:             ind_log++;
Chris@0:         }
Chris@0:         // fill negative values
Chris@0:         for (size_t ind = 0; ind < (m_warp_params.num_warps - 1) / 2; ind++) {
Chris@0:             m_warpings.chirp_rates[ind] = -m_warpings.chirp_rates[m_warp_params.num_warps - 1 - ind];
Chris@0:         }
Chris@0:     }
Chris@0: 
Chris@0:     // compute relative frequency deviation
Chris@0:     for (size_t i = 0; i < m_warpings.nsamps_torig; i++)
Chris@0:         for (size_t j = 0; j < m_warp_params.num_warps; j++)
Chris@0:             freq_relative[j * m_warpings.nsamps_torig + i] = 1.0 + t_orig[i] * m_warpings.chirp_rates[j];
Chris@0:     //freq_relative[i * m_warpings.nsamps_torig + j] = 1.0 + t_orig[i] * m_warpings.chirp_rates[j];
Chris@0:     //freq_relative[i][j] = 1.0 + t_orig[i] * m_warpings.chirp_rates[j];
Chris@0: 
Chris@0: 	printf("	End of define_warps_linear_chirps().\n");
Chris@0: 
Chris@0: }
Chris@0: 
Chris@0: void
Chris@0: FChTransformF0gram::design_LPF() {
Chris@0: 
Chris@0: 	printf("	Running design_LPF().\n");
Chris@0:     //    in = (fftw_complex*) fftw_malloc(sizeof (fftw_complex) * tamanoVentana);
Chris@0:     //    out = (fftw_complex*) fftw_malloc(sizeof (fftw_complex) * tamanoVentana);
Chris@0:     //    in_window = (float*) fftw_malloc(sizeof (float) * tamanoVentana);
Chris@0:     //    p = fftw_plan_dft_1d(tamanoVentana, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
Chris@0:     double *lp_LPFWindow_aux = new double[m_blockSize/2+1];
Chris@0:     mp_LPFWindow = new double[m_blockSize/2+1];
Chris@0:     
Chris@0:     size_t i_max = (size_t) ((2.0*m_fmax/m_fs) * ( (double)m_blockSize / 2.0 + 1.0 ));
Chris@0: 	//printf("	i_max = %d.\n", (int)i_max);
Chris@0:     for (size_t i = 0; i < m_blockSize/2+1; i++) {
Chris@0:         if (i >= i_max) {
Chris@0:             lp_LPFWindow_aux[i] = 0.0;
Chris@0:         } else {
Chris@0:             lp_LPFWindow_aux[i] = 1.0;
Chris@0:         }
Chris@0:     }
Chris@0:     LPF_time = (double*)fftw_malloc(sizeof ( double) * m_warpings.nsamps_torig);
Chris@0: 		//memset((char*)LPF_time, 0, m_warpings.nsamps_torig * sizeof(double));	
Chris@0: 			// sustituyo el memset por un for:
Chris@0: 			for (size_t i = 0; i < m_warpings.nsamps_torig; i++) {
Chris@0: 				LPF_time[i] = 0.0;
Chris@0: 			}
Chris@0: 		#ifdef DEBUG
Chris@0: 			printf("	Corrio primer memset...\n");
Chris@0: 		#endif
Chris@0:     LPF_frequency = (fftw_complex*)fftw_malloc(sizeof ( fftw_complex) * (m_warpings.nsamps_torig/2 + 1)); //tamaño de la fft cuando la entrada es real
Chris@0: 		//memset((char*)LPF_frequency, 0, sizeof(fftw_complex) * (m_warpings.nsamps_torig/2 + 1));
Chris@0: 			// sustituyo el memset por un for:
Chris@0: 			for (size_t i = 0; i < (m_warpings.nsamps_torig/2 + 1); i++) {
Chris@0: 				LPF_frequency[i][0] = 0.0;
Chris@0: 				LPF_frequency[i][1] = 0.0;
Chris@0: 			}
Chris@0: //	for (int i=0; i<(m_blockSize/2+1); i++) {
Chris@0: //		LPF_frequency[i] =  new fftw_complex;
Chris@0: //	}
Chris@0:     plan_forward_LPF = fftw_plan_dft_r2c_1d(m_blockSize, LPF_time, LPF_frequency, FFTW_ESTIMATE);
Chris@0:     plan_backward_LPF = fftw_plan_dft_c2r_1d(m_warpings.nsamps_torig, LPF_frequency, LPF_time, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
Chris@0:     
Chris@0: 	size_t winWidth = 11;
Chris@0:     double *lp_hanningWindow = new double[winWidth]; 
Chris@0:     double accum=0;
Chris@0:     for (size_t i = 0; i < winWidth; i++) {
Chris@0:         lp_hanningWindow[i]=0.5*(1.0-cos(2*M_PI*(double)(i+1)/((double)winWidth+1.0)));
Chris@0:         accum+=lp_hanningWindow[i];
Chris@0:         
Chris@0:     }
Chris@0:     for (size_t i = 0; i < winWidth; i++) { //window normalization
Chris@0:         lp_hanningWindow[i]=lp_hanningWindow[i]/accum;
Chris@0: 			//printf("	lp_hanningWindow[%d] = %f.\n",(int)i,lp_hanningWindow[i]);
Chris@0:     }
Chris@0:     for (size_t i = 0; i < m_blockSize/2+1; i++) {
Chris@0:         //if (((i-(winWidth-1)/2)<0)||(i+(winWidth-1))/2>m_blockSize/2-1) {//consideramos winWidth impar, si la ventana sale del arreglo se rellena con el valor origianl
Chris@0: 		if ( (i > (i_max + (winWidth-1)/2)) ||  (i <= (i_max - (winWidth-1)/2)) ) {
Chris@0:             mp_LPFWindow[i]=lp_LPFWindow_aux[i];
Chris@0: 			//printf("	entro al if en i=%d.\n",(int)i);
Chris@0:         } else {
Chris@0:             accum=0;
Chris@0:         	for (size_t j = -((winWidth-1)/2); j <= (winWidth-1)/2; j++) {
Chris@0:             	accum+=lp_LPFWindow_aux[i-j]*lp_hanningWindow[j+(winWidth-1)/2];
Chris@0: 				//printf("	accum = %f.\n",accum);
Chris@0:         	}
Chris@0:             mp_LPFWindow[i]=accum;
Chris@0:         }
Chris@0:     }
Chris@0: 
Chris@0:     delete[] lp_LPFWindow_aux;
Chris@0:     delete[] lp_hanningWindow;
Chris@0: 
Chris@0: 	printf("	End of design_LPF().\n");
Chris@0: }
Chris@0: 
Chris@0: void FChTransformF0gram::apply_LPF() {
Chris@0:     fftw_execute(plan_forward_LPF);
Chris@0:     for (size_t i = 0; i < m_blockSize/2+1; i++) {
Chris@0:         LPF_frequency[i][0]*=mp_LPFWindow[i];
Chris@0:         LPF_frequency[i][1]*=mp_LPFWindow[i];
Chris@0:     }
Chris@0:     fftw_execute(plan_backward_LPF);
Chris@0: 
Chris@0: 	// TODO ver si hay que hacer fftshift para corregir la fase respecto al centro del frame.
Chris@0: 	// nota: además de aplicar el LPF, esta función resamplea la señal original.
Chris@0: }
Chris@0: 
Chris@0: void FChTransformF0gram::clean_LPF() {
Chris@0:     delete[] mp_LPFWindow;
Chris@0: 
Chris@0: 	fftw_destroy_plan(plan_forward_LPF);
Chris@0: 	fftw_destroy_plan(plan_backward_LPF);
Chris@0: 	fftw_free(LPF_time);
Chris@0: 	fftw_free(LPF_frequency);
Chris@0: }
Chris@0: 
Chris@0: void FChTransformF0gram::reset() {
Chris@0: 
Chris@0:     printf("FUNCTION CALL: %s\n", __FUNCTION__);
Chris@0:     // Clear buffers, reset stored values, etc
Chris@0: 
Chris@0: 	delete [] m_warpings.pos_int;
Chris@0: 	delete [] m_warpings.pos_frac;
Chris@0: 
Chris@0:     fftw_destroy_plan(planFFT);
Chris@0: 	fftw_free(in); 
Chris@0: 	fftw_free(out);
Chris@0: 
Chris@0: 	clean_LPF();
Chris@0: 
Chris@0: 	delete [] m_timeWindow;
Chris@0: 
Chris@0: 	delete [] mp_HanningWindow;
Chris@0: 
Chris@0: 	// Warping
Chris@0: 	delete [] x_warping;
Chris@0: 	fftw_destroy_plan(plan_forward_xwarping);
Chris@0: 	fftw_free(m_absFanChirpTransform); 
Chris@0: 	fftw_free(m_auxFanChirpTransform);
Chris@0: 
Chris@0: 	// design_GLogS
Chris@0: 	delete [] m_glogs_f0;
Chris@0: 	delete [] m_glogs;
Chris@0: 	delete [] m_glogs_n;
Chris@0: 	delete [] m_glogs_index;
Chris@0: 	delete [] m_glogs_posint;
Chris@0: 	delete [] m_glogs_posfrac;
Chris@0: 	delete [] m_glogs_third_harmonic_posint;
Chris@0: 	delete [] m_glogs_third_harmonic_posfrac;
Chris@0: 	delete [] m_glogs_third_harmonic;
Chris@0: 	delete [] m_glogs_fifth_harmonic_posint;
Chris@0: 	delete [] m_glogs_fifth_harmonic_posfrac;
Chris@0: 	delete [] m_glogs_fifth_harmonic;
Chris@0: 	delete [] m_glogs_f0_preference_weights;
Chris@0: 	delete [] m_glogs_median_correction;
Chris@0: 	delete [] m_glogs_sigma_correction;
Chris@0: 	delete [] m_glogs_hf_smoothing_window;
Chris@0: 
Chris@0: }
Chris@0: 
Chris@0: FChTransformF0gram::FeatureSet
Chris@0: FChTransformF0gram::process(const float *const *inputBuffers, Vamp::RealTime timestamp) {
Chris@0: 
Chris@0:     printf("FUNCTION CALL: %s\n", __FUNCTION__);
Chris@0: 
Chris@0:     //    // Do actual work!
Chris@0:     //
Chris@0: 
Chris@0: 	/* PSEUDOCÓDIGO:
Chris@0: 	- Aplicar FFT al frame entero.	
Chris@0: 	- Filtro pasabajos en frecuencia.
Chris@0: 	- FFT inversa al frame entero.
Chris@0: -----------------------------------------------------------------------------
Chris@0: 	- Para cada warp: *Si es un espectrograma direccional (un solo warp 
Chris@0: 	  => no es para cada warp sino para el elegido)
Chris@0: 		- Hacer la interpolación con interp1q.
Chris@0: 		- Aplicar la FFT al frame warpeado.
Chris@0: 		- (Opcional) GLogS.
Chris@0: 	- ...
Chris@0: 	*/
Chris@0: 
Chris@0: //---------------------------------------------------------------------------
Chris@0: 	FeatureSet fs;
Chris@0: 	
Chris@0: 	#ifdef DEBUG
Chris@0: 		printf("\n	----- DEBUG INFORMATION ----- \n");
Chris@0: 		printf("	m_fs = %f Hz.\n",m_fs);
Chris@0: 		printf("	fs_orig = %f Hz.\n",m_warpings.fs_orig);
Chris@0: 		printf("	fs_warp = %f Hz.\n",m_warpings.fs_warp);
Chris@0: 		printf("	m_nfft = %d.\n",m_nfft);
Chris@0: 		printf("	m_blockSize = %d.\n",m_blockSize);
Chris@0: 		printf("	m_warpings.nsamps_torig = %d.\n",m_warpings.nsamps_torig);
Chris@0: 		printf("	m_warp_params.num_warps = %d.\n",m_warp_params.num_warps);
Chris@0: 		printf("	m_glogs_harmonic_count = %d.\n",m_glogs_harmonic_count);
Chris@0: 	#endif
Chris@0: 
Chris@0:     // size_t n = m_nfft/2 + 1;
Chris@0: 	// double *tbuf = in_window;
Chris@0: 
Chris@0: 	for (size_t i = 0; i < m_blockSize; i++) {
Chris@0:         LPF_time[i] = (double)(inputBuffers[0][i]) * m_timeWindow[i];
Chris@0:     }
Chris@0: 
Chris@0: //	#ifdef DEBUG
Chris@0: //		printf("	HASTA ACÁ ANDA!!!\n");
Chris@0: //		cout << flush;
Chris@0: //	#endif
Chris@0: 
Chris@0: 	apply_LPF();
Chris@0: 	// Señal filtrada queda en LPF_time
Chris@0: 
Chris@0: 	Feature feature;
Chris@0:     feature.hasTimestamp = false;
Chris@0: 
Chris@0: 
Chris@0: 	/* Solo a modo de prueba, voy a poner la salida del filtrado en «in» y
Chris@0: 	voy a mostrar la FFT de eso, para ver el efecto del filtrado. */
Chris@0: //    for (size_t i = 0; i < m_nfft; i++) {
Chris@0: //        in[i][0] = tbuf[i];
Chris@0: //        in[i][1] = 0;
Chris@0: //    }
Chris@0: //	fftw_execute(planFFT);
Chris@0: //	double real, imag;
Chris@0: //	for (size_t i=0; i<n; ++i) {		// preincremento?? ver version de nacho
Chris@0: //		real = out[i][0];
Chris@0: //		imag = out[i][1];
Chris@0: //		feature.values.push_back(real*real + imag*imag);
Chris@0: //	}
Chris@0: //	fs[0].push_back(feature);
Chris@0: 
Chris@0: //	float real; 
Chris@0: //	float imag;
Chris@0: //	for (size_t i=0; i<m_blockSize/2+1; i++) {
Chris@0: //		real = (float)(LPF_frequency[i][0]);
Chris@0: //		imag = (float)(LPF_frequency[i][1]);
Chris@0: //		feature.values.push_back(real*real+imag*imag);
Chris@0: //		//feature.values.push_back((float)(mp_LPFWindow[i]));
Chris@0: //	}
Chris@0: 
Chris@0: // ----------------------------------------------------------------------------------------------
Chris@0: // 		Hanning window & FFT for all warp directions
Chris@0: 
Chris@0: 	double max_glogs = -DBL_MAX;
Chris@0: 	size_t ind_max_glogs = 0;
Chris@0: 
Chris@0: 	for (size_t i_warp = 0; i_warp < m_warp_params.num_warps; i_warp++) {
Chris@0: 		// Interpolate
Chris@0: 		interp1q(LPF_time, (m_warpings.pos_int) + i_warp*m_warp_params.nsamps_twarp, m_warpings.pos_frac + i_warp*m_warp_params.nsamps_twarp, x_warping, m_warp_params.nsamps_twarp);
Chris@0: 
Chris@0: 		// Apply window
Chris@0: 		for (size_t i = 0; i < m_warp_params.nsamps_twarp; i++) {
Chris@0: 			x_warping[i] *= mp_HanningWindow[i];
Chris@0: 		}
Chris@0: 
Chris@0: 		// Transform
Chris@0: 		fftw_execute(plan_forward_xwarping);
Chris@0: 
Chris@0: 		// Copy result
Chris@0: 		//memcpy(m_absFanChirpTransform + i_warp*(m_warp_params.nsamps_twarp/2+1), m_auxFanChirpTransform, (m_warp_params.nsamps_twarp/2+1)*sizeof(fftw_complex)); asi como esta no funciona
Chris@0: 		double *aux_abs_fcht = m_absFanChirpTransform + i_warp*(m_warp_params.nsamps_twarp/2+1);
Chris@0: 		for (size_t i = 0; i < (m_warp_params.nsamps_twarp/2+1); i++) {
Chris@0: 			aux_abs_fcht[i] = log10(1.0 + 10.0*sqrt(m_auxFanChirpTransform[i][0]*m_auxFanChirpTransform[i][0]+m_auxFanChirpTransform[i][1]*m_auxFanChirpTransform[i][1]));
Chris@0: 			// smoothing high frequency values
Chris@0: 			//aux_abs_fcht[i] *= m_glogs_hf_smoothing_window[i];
Chris@0: 		}
Chris@0: 		
Chris@0: //      -----------------------------------------------------------------------------------------
Chris@0: // 		GLogS
Chris@0: 		interp1q(aux_abs_fcht, m_glogs_posint, m_glogs_posfrac, m_glogs_interp, m_glogs_harmonic_count);
Chris@0: 		size_t glogs_ind = 0;
Chris@0: 		for (size_t i = 0; i < m_glogs_num_f0s; i++) {
Chris@0: 			double glogs_accum = 0;
Chris@0: 			for (size_t j = 1; j <= m_glogs_n[i]; j++) {
Chris@0: 				glogs_accum += m_glogs_interp[glogs_ind++];
Chris@0: 			}
Chris@0: 			m_glogs[i + i_warp*m_glogs_num_f0s] = glogs_accum/(double)m_glogs_n[i];
Chris@0: 		}
Chris@0: 		//printf("	glogs_ind = %d.\n",(int)glogs_ind);
Chris@0: 
Chris@0: //		Sub/super harmonic correction
Chris@0: 		interp1q(m_glogs + i_warp*m_glogs_num_f0s, m_glogs_third_harmonic_posint, m_glogs_third_harmonic_posfrac, m_glogs_third_harmonic, (m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct);
Chris@0: 		interp1q(m_glogs + i_warp*m_glogs_num_f0s, m_glogs_fifth_harmonic_posint, m_glogs_fifth_harmonic_posfrac, m_glogs_fifth_harmonic, (m_f0_params.num_octs+1)*m_f0_params.num_f0s_per_oct);
Chris@0: 		for (size_t i = m_glogs_num_f0s-1; i >= m_glogs_init_f0s; i--) {
Chris@0: 			m_glogs[i + i_warp*m_glogs_num_f0s] -= MAX(MAX(m_glogs[i-m_f0_params.num_f0s_per_oct + i_warp*m_glogs_num_f0s],m_glogs_third_harmonic[i-m_glogs_init_f0s]),m_glogs_fifth_harmonic[i-m_glogs_init_f0s]);
Chris@0: 			//m_glogs[i] -= MAX(m_glogs[i-m_f0_params.num_f0s_per_oct],m_glogs_third_harmonic[i-m_glogs_init_f0s]);
Chris@0: 		}
Chris@0: 		for (size_t i = m_glogs_init_f0s; i < m_glogs_num_f0s-m_f0_params.num_f0s_per_oct; i++) {
Chris@0: 			m_glogs[i + i_warp*m_glogs_num_f0s] -= 0.3*m_glogs[i+m_f0_params.num_f0s_per_oct + i_warp*m_glogs_num_f0s];
Chris@0: 			// Median, sigma $ weights correction
Chris@0: 			m_glogs[i + i_warp*m_glogs_num_f0s] = (m_glogs[i + i_warp*m_glogs_num_f0s]-m_glogs_median_correction[i-m_glogs_init_f0s])*m_glogs_sigma_correction[i-m_glogs_init_f0s]*m_glogs_f0_preference_weights[i-m_glogs_init_f0s];
Chris@0: 		}
Chris@0: 	
Chris@0: 		// Look for maximum value to determine best direction
Chris@0: 		for (size_t i = m_glogs_init_f0s; i < m_glogs_num_f0s-m_f0_params.num_f0s_per_oct; i++) {
Chris@0: 			if (m_glogs[i + i_warp*m_glogs_num_f0s] > max_glogs) {
Chris@0: 				max_glogs = m_glogs[i + i_warp*m_glogs_num_f0s];
Chris@0: 				ind_max_glogs = i_warp;
Chris@0: 			}
Chris@0: 		}
Chris@0: 	}	
Chris@0: 	
Chris@0: // ----------------------------------------------------------------------------------------------
Chris@0: 
Chris@0: 	for (size_t i=m_glogs_init_f0s; i< m_glogs_num_f0s - m_f0_params.num_f0s_per_oct; i++) {
Chris@0: 	//for (size_t i=0; i<(m_warp_params.nsamps_twarp/2+1); i++) {
Chris@0: 		//feature.values.push_back((float)(m_warpings.pos_int[i])+ (float)(m_warpings.pos_frac[i]));
Chris@0: 		//feature.values.push_back((float)(phi[i]*100000.0));
Chris@0: 		//feature.values.push_back((float)(t_orig[i]));
Chris@0: 		//feature.values.push_back((float)(pos1[i]));
Chris@0: 		//feature.values.push_back((float)x_warping[i]);
Chris@0: 		//feature.values.push_back(m_absFanChirpTransform[i + ind_max_glogs*(m_warp_params.nsamps_twarp/2+1)]);
Chris@0: 		//feature.values.push_back((float)m_glogs[i+(long)ind_max_glogs*(long)m_glogs_num_f0s]);
Chris@0: 		switch (m_f0gram_mode) {
Chris@0: 			case 1:		
Chris@0: 					max_glogs = -DBL_MAX;
Chris@0: 					for	(size_t i_warp = 0; i_warp < m_warp_params.num_warps; i_warp++) {
Chris@0: 						if (m_glogs[i + i_warp*m_glogs_num_f0s] > max_glogs) {
Chris@0: 							max_glogs = m_glogs[i + i_warp*m_glogs_num_f0s];
Chris@0: 							ind_max_glogs = i_warp;
Chris@0: 						}
Chris@0: 					}
Chris@0: 					feature.values.push_back((float)max_glogs);
Chris@0: 					break;
Chris@0: 			case 0:
Chris@0: 					feature.values.push_back((float)m_glogs[i+(size_t)ind_max_glogs*(size_t)m_glogs_num_f0s]);
Chris@0: 					break;
Chris@0: 		}
Chris@0: 		//feature.values.push_back((float)m_glogs_hf_smoothing_window[i]);
Chris@0: 	}
Chris@0: 
Chris@0: // ----------------------------------------------------------------------------------------------
Chris@0: 
Chris@0: 	fs[0].push_back(feature);
Chris@0: 
Chris@0: 	#ifdef DEBUG
Chris@0: 		printf("	----------------------------- \n");
Chris@0: 	#endif
Chris@0: 	
Chris@0: 	return fs;
Chris@0: //---------------------------------------------------------------------------
Chris@0: 
Chris@0:     //return FeatureSet();
Chris@0: }
Chris@0: 
Chris@0: FChTransformF0gram::FeatureSet
Chris@0: FChTransformF0gram::getRemainingFeatures() {
Chris@0: 
Chris@0:     printf("FUNCTION CALL: %s\n", __FUNCTION__);
Chris@0: 
Chris@0:     return FeatureSet();
Chris@0: }
Chris@0: 
Chris@0: void
Chris@0: FChTransformF0gram::design_time_window() {
Chris@0: 
Chris@0: 	printf("	Running design_time_window().\n");
Chris@0: 
Chris@0: 	size_t transitionWidth = (size_t)m_blockSize/128 + 1;;
Chris@0:     m_timeWindow = new double[m_blockSize];
Chris@0: 	double *lp_transitionWindow = new double[transitionWidth];
Chris@0: 
Chris@0: 	//memset(m_timeWindow, 1.0, m_blockSize);
Chris@0: 	for (size_t i = 0; i < m_blockSize; i++) {
Chris@0: 		m_timeWindow[i] = 1.0;
Chris@0: 	}
Chris@0: 
Chris@0: 	for (size_t i = 0; i < transitionWidth; i++) {
Chris@0:         lp_transitionWindow[i]=0.5*(1.0-cos(2*M_PI*(double)(i+1)/((double)transitionWidth+1.0)));
Chris@0:     }
Chris@0: 
Chris@0: 	for (size_t i = 0; i < transitionWidth/2; i++) {
Chris@0: 		m_timeWindow[i] = lp_transitionWindow[i];
Chris@0: 		m_timeWindow[m_blockSize-1-i] = lp_transitionWindow[transitionWidth-1-i];
Chris@0: 	}
Chris@0: 
Chris@0: 	#ifdef DEBUG
Chris@0: 		for (int i = 0; i < m_blockSize; i++) {
Chris@0: 			if ((i<transitionWidth)) {
Chris@0: 				printf("	m_timeWindow[%d] = %f.\n",i,m_timeWindow[i]);
Chris@0: 			}
Chris@0: 		}
Chris@0: 	#endif
Chris@0: 
Chris@0: 	delete [] lp_transitionWindow;
Chris@0: 
Chris@0: 	printf("	End of design_time_window().\n");
Chris@0: }
Chris@0: