mathieu@0: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
mathieu@0: 
mathieu@0: /*
mathieu@0:  Calcium Signal Analyser Vamp Plugin
mathieu@0:  
mathieu@2:  Transient detection and characterisation for calcium signals.
mathieu@0:  
mathieu@2:  Centre for Digital Music, Queen Mary University of London.
mathieu@2:  This file copyright 2010-2011 Mathieu Barthet and QMUL.
mathieu@0:  
mathieu@2:  Based on the QM Vamp note onset detector plugin by Christian 
mathieu@2:  Landone, Chris Duxbury, and Juan Pablo Bello.
mathieu@0:  
mathieu@2:  This program is free software; you can redistribute it and/or
mathieu@2:  modify it under the terms of the GNU General Public License as
mathieu@2:  published by the Free Software Foundation; either version 2 of the
mathieu@2:  License, or (at your option) any later version.  See the file
mathieu@2:  COPYING included with this distribution for more information.
mathieu@2:   
mathieu@2:  Version: 2
mathieu@2: */
mathieu@0: 
mathieu@0: #include "CalciumSignalAnalyser.h"
mathieu@0: 
mathieu@0: #include <dsp/onsets/PeakPicking.h>
mathieu@0: #include <dsp/tempotracking/TempoTrack.h>
mathieu@0: 
mathieu@0: using std::string;
mathieu@0: using std::vector;
mathieu@0: using std::cerr;
mathieu@0: using std::endl;
mathieu@0: 
mathieu@0: #define MAX_PEAKFREQ 10; //maximum peak frequency in Hz
mathieu@0: 
mathieu@0: //#define VERBOSE //remove comment to have output messages for data display or debugging information
mathieu@0: 
mathieu@0: CalciumSignalAnalyser::CalciumSignalAnalyser(float inputSampleRate) :
mathieu@0: Vamp::Plugin(inputSampleRate),
mathieu@0: m_inputSampleRate(inputSampleRate),
mathieu@0: m_blockSize(0),
mathieu@0: m_stepSize(0),
mathieu@0: data(0),
mathieu@0: time(0),
mathieu@0: m_sensitivity(50),
mathieu@0: m_delta(0.2),
mathieu@2: m_mfwindur(10),
mathieu@0: processcounter(0)
mathieu@0: {
mathieu@0: #ifdef VERBOSE
mathieu@0: 	cout << "Sample rate: " << inputSampleRate << endl;
mathieu@0: #endif
mathieu@0: }
mathieu@0: 
mathieu@0: CalciumSignalAnalyser::~CalciumSignalAnalyser()
mathieu@0: {	
mathieu@0: #ifdef VERBOSE
mathieu@0: 	cout << "Destructor" << endl;
mathieu@0: #endif
mathieu@0: 	
mathieu@0: 	if (!data.empty()){
mathieu@0: #ifdef VERBOSE
mathieu@0: 		cout << "Data vector reset" << endl;
mathieu@0: #endif
mathieu@0: 		data.clear(); //clears the vector's elements
mathieu@0: 	}
mathieu@0: 	
mathieu@0: 	if (!time.empty()){
mathieu@0: #ifdef VERBOSE
mathieu@0: 		cout << "Time vector reset" << endl;
mathieu@0: #endif
mathieu@0: 		time.clear(); //clears the vector's elements
mathieu@0: 	}
mathieu@0: }
mathieu@0: 
mathieu@0: string
mathieu@0: CalciumSignalAnalyser::getIdentifier() const
mathieu@0: {
mathieu@0:     return "peakdetection";
mathieu@0: }
mathieu@0: 
mathieu@0: string
mathieu@0: CalciumSignalAnalyser::getName() const
mathieu@0: {
mathieu@0:     return "Calcium Signal Analyser";
mathieu@0: }
mathieu@0: 
mathieu@0: string
mathieu@0: CalciumSignalAnalyser::getDescription() const
mathieu@0: {
mathieu@0:     return "Automatic detection and characterisation of the signal's transients";
mathieu@0: }
mathieu@0: 
mathieu@0: string
mathieu@0: CalciumSignalAnalyser::getMaker() const
mathieu@0: {
mathieu@0:     return "Mathieu Barthet";
mathieu@0: }
mathieu@0: 
mathieu@0: int
mathieu@0: CalciumSignalAnalyser::getPluginVersion() const
mathieu@0: {
mathieu@2:     return 2;
mathieu@0: }
mathieu@0: 
mathieu@0: string
mathieu@0: CalciumSignalAnalyser::getCopyright() const
mathieu@0: {
mathieu@0:     return "Plugin by Mathieu Barthet. Based on the note onset detector plugin by Christian Landone, Chris Duxbury and Juan Pablo Bello. Copyright (c) 2010 QMUL - All Rights Reserved";
mathieu@0: }
mathieu@0: 
mathieu@0: CalciumSignalAnalyser::ParameterList
mathieu@0: CalciumSignalAnalyser::getParameterDescriptors() const
mathieu@0: {
mathieu@0:     ParameterList list;
mathieu@0: 	ParameterDescriptor desc;
mathieu@0: 	
mathieu@0:     desc.identifier = "sensitivity";
mathieu@0:     desc.name = "Peak-picking sensitivity";
mathieu@0:     desc.description = "Sensitivity of the peak detection (criterion based on quadratic interpolation)";
mathieu@0:     desc.minValue = 0;
mathieu@0:     desc.maxValue = 100;
mathieu@0:     desc.defaultValue = 50;
mathieu@0:     desc.isQuantized = true;
mathieu@0:     desc.quantizeStep = 1;
mathieu@0:     desc.unit = "%";
mathieu@0:     list.push_back(desc);
mathieu@0: 	
mathieu@0: 	desc.identifier = "deltathreshold";
mathieu@0:     desc.name = "Peak-picking offset threshold (delta)";
mathieu@0:     desc.description = "Offset threshold aiming at improving the peak detection by discarding the noise";
mathieu@0:     desc.minValue = 0;
mathieu@0:     desc.maxValue = 1;
mathieu@0:     desc.defaultValue = 0.2;
mathieu@0:     desc.isQuantized = true;
mathieu@0:     desc.quantizeStep = 0.05;
mathieu@0:     desc.unit = "";
mathieu@0:     list.push_back(desc);
mathieu@0: 	
mathieu@2: 	desc.identifier = "mfwindowduration";
mathieu@2:     desc.name = "Duration of the median filtering window";
mathieu@2:     desc.description = "Sets the duration of the median filtering window";
mathieu@2:     desc.minValue = 2;
mathieu@2:     desc.maxValue = 100;
mathieu@2:     desc.defaultValue = 10;
mathieu@2:     desc.isQuantized = true;
mathieu@2:     desc.quantizeStep = 0.1;
mathieu@2:     desc.unit = "s";
mathieu@2:     list.push_back(desc);
mathieu@2: 
mathieu@0:     return list;
mathieu@0: }
mathieu@0: 
mathieu@0: float
mathieu@0: CalciumSignalAnalyser::getParameter(std::string name) const
mathieu@0: {	
mathieu@0: 	if (name == "sensitivity"){
mathieu@0: 		return m_sensitivity;
mathieu@0: 	}
mathieu@0: 	else if (name == "deltathreshold") {
mathieu@0: 		return m_delta;	
mathieu@2: 	}
mathieu@2: 	else if (name == "mfwindowduration") {
mathieu@2: 		return m_mfwindur;	
mathieu@0: 	}	
mathieu@0: 	
mathieu@0:     return 0.0;
mathieu@0: }
mathieu@0: 
mathieu@0: void
mathieu@0: CalciumSignalAnalyser::setParameter(std::string name, float value)
mathieu@0: {	
mathieu@0: 	if (name == "sensitivity") {
mathieu@0:         if (m_sensitivity == value) return;
mathieu@0:         m_sensitivity = value;
mathieu@0:         
mathieu@0: 	}
mathieu@0: 	else if (name == "deltathreshold") {
mathieu@0: 		if (m_delta == value) return;
mathieu@0:         m_delta = value;
mathieu@0: 	}
mathieu@2: 	else if (name == "mfwindowduration") {
mathieu@2: 		if (m_mfwindur == value) return;
mathieu@2:         m_mfwindur = value;
mathieu@2: 	}
mathieu@0: }
mathieu@0: 
mathieu@0: 
mathieu@0: CalciumSignalAnalyser::ProgramList
mathieu@0: CalciumSignalAnalyser::getPrograms() const
mathieu@0: {
mathieu@0:     ProgramList programs;
mathieu@0: 
mathieu@0:     return programs;
mathieu@0: }
mathieu@0: 
mathieu@0: std::string
mathieu@0: CalciumSignalAnalyser::getCurrentProgram() const
mathieu@0: {	
mathieu@0: 	return "";
mathieu@0: }
mathieu@0: 
mathieu@0: void
mathieu@0: CalciumSignalAnalyser::selectProgram(std::string program)
mathieu@0: {
mathieu@0: }
mathieu@0: 
mathieu@0: bool
mathieu@0: CalciumSignalAnalyser::initialise(size_t channels, size_t stepSize, size_t blockSize)
mathieu@0: {
mathieu@0: 	
mathieu@0: #ifdef VERBOSE
mathieu@0: 	cout << "Initialise" << endl;
mathieu@0: #endif
mathieu@0: 
mathieu@0: 	m_blockSize = 1;
mathieu@0: 	m_stepSize = 1;
mathieu@0: 	
mathieu@0: #ifdef VERBOSE
mathieu@0: 	cout << m_blockSize << " " << m_stepSize << endl;
mathieu@0: #endif
mathieu@0: 	
mathieu@0:     if (channels < getMinChannelCount() ||
mathieu@0: 		channels > getMaxChannelCount()) {
mathieu@0:         std::cerr << "CalciumSignalAnalyser::initialise: Unsupported channel count: "
mathieu@0: 		<< channels << std::endl;
mathieu@0:         return false;
mathieu@0:     }
mathieu@0: 	
mathieu@0:     if (stepSize != getPreferredStepSize()) {
mathieu@0: 		
mathieu@0: 		std::cerr << "ERROR: CalciumSignalAnalyser::initialise: the step size has to be 1." << std::endl;
mathieu@0: 		return false;
mathieu@0: 	}
mathieu@0: 	
mathieu@0:     if (blockSize != getPreferredBlockSize()) {
mathieu@0: 		
mathieu@0: 		std::cerr << "ERROR: CalciumSignalAnalyser::initialise: the block size has to be 1." << std::endl;
mathieu@0: 		return false;
mathieu@0: 	}
mathieu@0: 	
mathieu@0:     return true;
mathieu@0: }
mathieu@0: 
mathieu@0: void
mathieu@0: CalciumSignalAnalyser::reset()
mathieu@0: {	
mathieu@0: #ifdef VERBOSE
mathieu@0: 	cout << "Reset" << endl;
mathieu@0: #endif
mathieu@0: 	
mathieu@0: 	if (!data.empty()){
mathieu@0: #ifdef VERBOSE
mathieu@0: 		cout << "Data vector reset" << endl;
mathieu@0: #endif
mathieu@0: 		data.clear(); //clears the vector's elements
mathieu@0: 	}
mathieu@0: 	
mathieu@0: 	if (!time.empty()){
mathieu@0: #ifdef VERBOSE
mathieu@0: 		cout << "Time vector reset" << endl;
mathieu@0: #endif
mathieu@0: 		time.clear(); //clears the vector's elements
mathieu@0: 	}
mathieu@0: }
mathieu@0: 
mathieu@0: size_t
mathieu@0: CalciumSignalAnalyser::getPreferredStepSize() const
mathieu@0: {	
mathieu@0: 	return getPreferredBlockSize();
mathieu@0: }
mathieu@0: 
mathieu@0: size_t
mathieu@0: CalciumSignalAnalyser::getPreferredBlockSize() const
mathieu@0: {
mathieu@0: 	return 1;
mathieu@0: }
mathieu@0: 
mathieu@0: size_t
mathieu@0: CalciumSignalAnalyser::getMinChannelCount() const
mathieu@0: {
mathieu@0:     return 1;
mathieu@0: }
mathieu@0: 
mathieu@0: size_t
mathieu@0: CalciumSignalAnalyser::getMaxChannelCount() const
mathieu@0: {
mathieu@0:     return 1;
mathieu@0: }
mathieu@0: 
mathieu@0: 
mathieu@0: CalciumSignalAnalyser::OutputList
mathieu@0: CalciumSignalAnalyser::getOutputDescriptors() const
mathieu@0: {	
mathieu@0: #ifdef VERBOSE
mathieu@0: 	cout << "get output desc" << endl;
mathieu@0: #endif
mathieu@0: 	
mathieu@0:     OutputList list;
mathieu@0:   	
mathieu@0:     OutputDescriptor pt;
mathieu@0:     pt.identifier = "peaktimes";
mathieu@0:     pt.name = "Peaks times";
mathieu@0:     pt.description = "Time positions of the peaks (in seconds)";
mathieu@0:     pt.unit = "";
mathieu@0:     pt.hasFixedBinCount = true;
mathieu@0:     pt.binCount = 0; //0 for time only feature
mathieu@0:     pt.sampleType = OutputDescriptor::VariableSampleRate;
mathieu@0: 	pt.hasDuration = false;
mathieu@0: 	pt.hasKnownExtents = false;
mathieu@0: 	pt.isQuantized = false;
mathieu@0:     pt.sampleRate = 0; //no sensible resolution, the time of the feature are given in the timestamp of the feature
mathieu@0: 	list.push_back(pt);
mathieu@0: 	
mathieu@0: 	OutputDescriptor po;
mathieu@0:     po.identifier = "peakonsets";
mathieu@0:     po.name = "Peaks onsets";
mathieu@0:     po.description = "Time positions of the peak onsets (in seconds)";
mathieu@0:     po.unit = "";
mathieu@0:     po.hasFixedBinCount = true;
mathieu@0:     po.binCount = 0; //0 for time only feature
mathieu@0:     po.sampleType = OutputDescriptor::VariableSampleRate;
mathieu@0: 	po.hasDuration = false;
mathieu@0: 	po.hasKnownExtents = false;
mathieu@0: 	po.isQuantized = false;
mathieu@0:     po.sampleRate = 0; //no sensible resolution, the time of the feature are given in the timestamp of the feature
mathieu@0: 	list.push_back(po);
mathieu@0: 	
mathieu@0: 	OutputDescriptor sdf;
mathieu@0:     sdf.identifier = "smoothed_df";
mathieu@0:     sdf.name = "Smoothed Detection Function";
mathieu@0:     sdf.description = "Smoothed probability function used for peak-picking";
mathieu@0:     sdf.unit = "";
mathieu@0:     sdf.hasFixedBinCount = true;
mathieu@0:     sdf.binCount = 1;
mathieu@0:     sdf.hasKnownExtents = false;
mathieu@0:     sdf.isQuantized = false;
mathieu@0:     sdf.sampleType = OutputDescriptor::FixedSampleRate;
mathieu@0: 	sdf.sampleRate = m_inputSampleRate; //the smoothed detection function has the same sample rate as the input signal
mathieu@0:     list.push_back(sdf);
mathieu@0: 		
mathieu@0: 	OutputDescriptor pf;
mathieu@0:     pf.identifier = "peak_frequency";
mathieu@0:     pf.name = "Peak frequency";
mathieu@0:     pf.description = "Frequency of the peaks (in Hz)";
mathieu@0:     //pf.unit = "Hz"; //if unit is Hz then the feature is treated as a pitch which is not relevant in this case
mathieu@0:     pf.unit = "";
mathieu@0: 	pf.hasFixedBinCount = true;
mathieu@0:     pf.binCount = 1; //1 corresponds to time and value
mathieu@0:     pf.sampleType = OutputDescriptor::VariableSampleRate;
mathieu@0: 	pf.sampleRate = 0; //no sensible resolution, the time of the feature is given in the timestamp of the feature
mathieu@0: 	pf.hasKnownExtents = true;
mathieu@0: 	pf.minValue = 0;
mathieu@0: 	pf.maxValue = MAX_PEAKFREQ;	
mathieu@0: 	pf.isQuantized = false;
mathieu@0: 	pf.hasDuration = true; //duration of the signal
mathieu@0: 	list.push_back(pf);
mathieu@0:     
mathieu@0: 	OutputDescriptor mipi;
mathieu@0:     mipi.identifier = "mean_interpeakinterval";
mathieu@0:     mipi.name = "Mean Inter Peak Interval";
mathieu@0:     mipi.description = "Mean of the time intervals between peaks (in seconds)";
mathieu@0:     mipi.unit = "s";
mathieu@0: 	mipi.hasFixedBinCount = true;
mathieu@0:     mipi.binCount = 1; //1 corresponds to time and value
mathieu@0:     mipi.sampleType = OutputDescriptor::VariableSampleRate;
mathieu@0: 	mipi.sampleRate = 0; //no sensible resolution, the time of the feature is given in the timestamp of the feature
mathieu@0: 	mipi.hasKnownExtents = false;
mathieu@0: 	mipi.isQuantized = false;
mathieu@0: 	mipi.hasDuration = true; //duration of the signal
mathieu@0: 	list.push_back(mipi);
mathieu@0: 	
mathieu@0:     return list;
mathieu@0: }
mathieu@0: 
mathieu@0: CalciumSignalAnalyser::FeatureSet
mathieu@0: CalciumSignalAnalyser::process(const float *const *inputBuffers,
mathieu@0:                        Vamp::RealTime timestamp)
mathieu@0: {	
mathieu@0: 	FeatureSet fs;
mathieu@0: 	
mathieu@0: 	if (m_blockSize == 0 || m_stepSize == 0) {
mathieu@0: 		cerr << "ERROR: CalciumSignalAnalyser::process: "
mathieu@0: 		<< "CalciumSignalAnalyser has not been initialised."
mathieu@0: 		<< endl;
mathieu@0: 		return fs;
mathieu@0:     }
mathieu@0: 	
mathieu@0: 	data.push_back(inputBuffers[0][0]); //assumes that only one sample is read at a time
mathieu@0: 	time.push_back(timestamp);
mathieu@0: 	
mathieu@0: 	processcounter ++;
mathieu@0: 	
mathieu@0: 	fs = FeatureSet(); //empty feature set
mathieu@0: 	
mathieu@0: 	return fs;	
mathieu@0: }
mathieu@0: 
mathieu@0: CalciumSignalAnalyser::FeatureSet
mathieu@0: CalciumSignalAnalyser::getRemainingFeatures()
mathieu@0: {
mathieu@0:    	FeatureSet returnFeatures;
mathieu@0: 	returnFeatures = FeatureSet(); 
mathieu@0: 	
mathieu@0: #ifdef VERBOSE
mathieu@0: 	cout << "Get remaining features" << endl << endl;
mathieu@0: #endif
mathieu@0: 	
mathieu@0: 	/***************************************
mathieu@0: 	 Peak times
mathieu@0: 	 ****************************************/
mathieu@0:     
mathieu@0: 	double aCoeffs[] = { 1.0000, -0.5949, 0.2348 };
mathieu@0:     double bCoeffs[] = { 0.1600,  0.3200, 0.1600 };
mathieu@0: 	
mathieu@0:     //The signal is directly treated as a detection function for the peak picking stage
mathieu@0: 	
mathieu@0:     PPickParams ppParams;
mathieu@0:     ppParams.length = data.size(); //length of the signal
mathieu@0: 	
mathieu@0: #ifdef VERBOSE
mathieu@0: 	cout << "Length: " << ppParams.length << endl;
mathieu@0: #endif
mathieu@0: 	
mathieu@0: 	// tau and cutoff appear to be unused in PeakPicking, but I've
mathieu@0:     // inserted some moderately plausible values rather than leave
mathieu@0:     // them unset.  The QuadThresh values come from trial and error.
mathieu@0:     // The rest of these are copied from ttParams in the BeatTracker
mathieu@0:     // code: I don't claim to know whether they're good or not --cc
mathieu@0:         
mathieu@0: 	ppParams.tau = m_stepSize / m_inputSampleRate; //not used for peak picking 
mathieu@0: 	
mathieu@0: 	ppParams.alpha = 9;
mathieu@0:     ppParams.cutoff = m_inputSampleRate/4; //not used for peak picking 
mathieu@0:     ppParams.LPOrd = 2;
mathieu@0:     ppParams.LPACoeffs = aCoeffs;
mathieu@0:     ppParams.LPBCoeffs = bCoeffs;
mathieu@0:     
mathieu@2: 	int NWin;
mathieu@2: 	float samples = m_inputSampleRate*m_mfwindur;
mathieu@2: 	if (samples < 4){
mathieu@2: 		NWin = 4; //minimal possible length of the window in samples
mathieu@2: #ifdef VERBOSE
mathieu@2: 		cout << "WARNING: CalciumSignalAnalyser::getRemainingFeatures:"
mathieu@2: 		<< "The duration of the median filtering window is too small considering the"
mathieu@2: 		<< "sample frequency of the signal. Duration of "<< NWin/m_inputSampleRate
mathieu@2: 		<< " s used instead."
mathieu@2: 		<< endl;
mathieu@2: #endif
mathieu@2: 	}else{
mathieu@2: 		NWin = floor(m_inputSampleRate*m_mfwindur);
mathieu@2: 		if (NWin%2 != 0) NWin--;
mathieu@2: 	}	
mathieu@2: 		
mathieu@2: 	ppParams.WinT.post = NWin/2;
mathieu@2:     ppParams.WinT.pre = NWin/2 - 1;
mathieu@0:     		
mathieu@0: 	ppParams.QuadThresh.a = (100 - m_sensitivity) / 1000.0;
mathieu@0:     ppParams.QuadThresh.b = 0;
mathieu@0:     ppParams.QuadThresh.c = (100 - m_sensitivity) / 1500.0;
mathieu@0: 	ppParams.delta = m_delta; //delta threshold used as an offset when computing the smoothed detection function
mathieu@0: 			
mathieu@0:     PeakPicking peakPicker(ppParams);
mathieu@0: 	
mathieu@0: 	double *ppSrc = new double[ppParams.length];
mathieu@0: 	if (ppSrc==NULL)	{
mathieu@0: 		cerr << "ERROR: CalciumSignalAnalyser::getRemainingFeatures:"
mathieu@0: 		<< "Unable to allocate memory (ppSrc)."
mathieu@0: 		<< endl;
mathieu@0: 		return returnFeatures;
mathieu@0: 	}
mathieu@0: 		
mathieu@0:     for (unsigned int i = 0; i < ppParams.length; ++i) {
mathieu@0:         
mathieu@0: 		ppSrc[i] = (double) data[i];
mathieu@0: 	}	
mathieu@0: 	
mathieu@0:     vector<int> peaks;
mathieu@0: 	vector<int> peak_frame;
mathieu@0: 	vector<int> onsets_frame;
mathieu@0: 	
mathieu@0: 	peakPicker.process(ppSrc, ppParams.length, peaks);
mathieu@0: 	
mathieu@0: 	for (size_t i = 0; i < peaks.size(); ++i) {
mathieu@0: 		
mathieu@0: 		//peak times
mathieu@0:         size_t index = peaks[i];
mathieu@0: 		size_t peakframe = peaks[i] * m_stepSize;
mathieu@0: 		peak_frame.push_back(peakframe);
mathieu@0: 		
mathieu@0: 		Feature featurepeaks;
mathieu@0: 		featurepeaks.hasTimestamp = true;	
mathieu@0: 		featurepeaks.timestamp = Vamp::RealTime::frame2RealTime(peakframe, lrintf(m_inputSampleRate));	
mathieu@0: 		
mathieu@0: 		returnFeatures[0].push_back(featurepeaks); // peak times are output 0
mathieu@0: 		
mathieu@0: 		//algorithm to detect the peak onsets
mathieu@0: 		 
mathieu@0: 		double prevDiff = 0.0;
mathieu@0: 		while (index > 1) {
mathieu@0: 			double diff = ppSrc[index] - ppSrc[index-1];
mathieu@0: 			if (diff < prevDiff * 0.9) break;
mathieu@0: 			prevDiff = diff;
mathieu@0: 			--index;
mathieu@0: 		}	
mathieu@0: 			
mathieu@0: 		size_t onsetframe = index * m_stepSize;		
mathieu@0: 		onsets_frame.push_back(onsetframe);
mathieu@0: 		
mathieu@0: 		Feature featureonsets;
mathieu@0: 		featureonsets.hasTimestamp = true;	
mathieu@0: 		featureonsets.timestamp = Vamp::RealTime::frame2RealTime(onsetframe, lrintf(m_inputSampleRate));	
mathieu@0: 		
mathieu@0: #ifdef VERBOSE
mathieu@0: 		cout << featureonsets.timestamp << endl;
mathieu@0: #endif
mathieu@0: 		returnFeatures[1].push_back(featureonsets); // peak onsets are output 1
mathieu@0:     }
mathieu@0: 	
mathieu@0:     for (unsigned int i = 0; i < ppParams.length; ++i) {
mathieu@0:         
mathieu@0:         Feature feature;
mathieu@0: 		feature.hasTimestamp = true;
mathieu@0: 		size_t frame = i * m_stepSize;
mathieu@0: 				
mathieu@0: 		feature.timestamp = Vamp::RealTime::frame2RealTime(frame, lrintf(m_inputSampleRate));
mathieu@0: 		
mathieu@0:         feature.values.push_back(ppSrc[i]);
mathieu@0:      	
mathieu@0: 		returnFeatures[2].push_back(feature); // smoothed df is output 2
mathieu@0:     }
mathieu@0: 	
mathieu@0: 	/***************************************
mathieu@0: 	 Peak frequency
mathieu@0: 	 ****************************************/
mathieu@0: 	
mathieu@0: 	double peakfreq, duration;
mathieu@0: 	int npeaks; //total number of peaks
mathieu@0: 	
mathieu@0: 	//last element of the time vector is the duration of the signal - converts the RealTime to a float
mathieu@0: 	duration = time.back()/Vamp::RealTime(1,0); //division by 1 s to obtain a double
mathieu@0: 		
mathieu@0: 	npeaks = peaks.size();
mathieu@0: 	peakfreq = npeaks/duration;
mathieu@0: 	
mathieu@0: #ifdef VERBOSE
mathieu@0: 	cout << npeaks << endl;
mathieu@0: 	cout << duration << endl;	
mathieu@0: 	cout << peakfreq << endl;
mathieu@0: #endif
mathieu@0: 	
mathieu@0: 	Feature pf;
mathieu@0: 	pf.hasTimestamp = true;
mathieu@0: 	pf.timestamp = time[0];
mathieu@0: 	pf.hasDuration = true;			
mathieu@0: 	pf.duration = time.back()-time[0];
mathieu@0: 	
mathieu@0: #ifdef VERBOSE
mathieu@0: 	cout << "processcounter" << processcounter << endl;
mathieu@0: 	
mathieu@0: 	cout << "time.back " << time.back().toText() << endl;
mathieu@0: 	cout << "time.size " << time.size() << endl;
mathieu@0: 	cout << "data.size " << data.size() << endl;
mathieu@0: 	cout << "time[0, 1 ...] " << time[0].toText() << " " << time[1].toText() << " " << time[2].toText() << " " << time[time.size()].toText() << endl;
mathieu@0: 	
mathieu@0: 	cout << "time[999] " << time[999].toText() << endl;
mathieu@0: 	cout << "data[999] " << data[999] << endl;
mathieu@0: 	cout << "time[1000] " << time[1000].toText() << endl;
mathieu@0: 	cout << "data[1000] " << data[1000]<< endl;
mathieu@0: #endif
mathieu@0: 	
mathieu@0: 	pf.values.push_back(peakfreq);
mathieu@0: 		
mathieu@0: 	char pflabel[256];
mathieu@0: 	sprintf(pflabel,"%f Hz",peakfreq);
mathieu@0: 	pf.label = pflabel;
mathieu@0: 	
mathieu@0: #ifdef VERBOSE
mathieu@0: 	cout << "label: " << pflabel << endl;
mathieu@0: #endif
mathieu@0: 	
mathieu@0: 	returnFeatures[3].push_back(pf); //peak frequency is output 3
mathieu@0: 	
mathieu@0: 	/***************************************
mathieu@0: 	 Mean Inter Peak Interval
mathieu@0: 	 ****************************************/
mathieu@0: 	
mathieu@0: 	double o1,o2;
mathieu@0: 	
mathieu@0: 	double sumipi_s = 0;
mathieu@0: 	double meanipi_s;
mathieu@0: 	
mathieu@0: 	for (int j = 0; j < npeaks-1; ++j) {
mathieu@0: 		o2 = Vamp::RealTime::frame2RealTime(peak_frame[j+1], lrintf(m_inputSampleRate))/Vamp::RealTime(1,0);
mathieu@0: 		o1 = Vamp::RealTime::frame2RealTime(peak_frame[j], lrintf(m_inputSampleRate))/Vamp::RealTime(1,0);
mathieu@0: #ifdef VERBOSE
mathieu@0: 		cout << "o1: " << o1 << endl;
mathieu@0: 		cout << "o2: " << o2 << endl;
mathieu@0: #endif
mathieu@0: 		sumipi_s = sumipi_s + o2 - o1;	
mathieu@0: 	}
mathieu@0: 	
mathieu@0: 	meanipi_s = sumipi_s/(npeaks-1);
mathieu@0: 	
mathieu@0: #ifdef VERBOSE
mathieu@0: 	cout << "Sum IPI: " << sumipi_s << endl;
mathieu@0: 	cout << "Mean IPI: " << meanipi_s << endl;
mathieu@0: #endif
mathieu@0: 	
mathieu@0: 	Feature mipi;
mathieu@0: 	mipi.hasTimestamp = true;
mathieu@0: 	mipi.timestamp = time[0];
mathieu@0: 	mipi.hasDuration = true;			
mathieu@0: 	mipi.duration = time.back()-time[0];
mathieu@0: 	
mathieu@0: 	mipi.values.push_back(meanipi_s); //mean inter peak interval in secs
mathieu@0: 	
mathieu@0: 	returnFeatures[4].push_back(mipi); //mean inter peak interval is output 4
mathieu@0: 	
mathieu@0: 	//clear data
mathieu@0: 	
mathieu@0: 	delete [] ppSrc;
mathieu@0: 	
mathieu@0: 	if (!data.empty()){
mathieu@0: #ifdef VERBOSE
mathieu@0: 		cout << "Data vector reset" << endl;
mathieu@0: #endif
mathieu@0: 		data.clear(); //clears the vector's elements
mathieu@0: 	}
mathieu@0: 	
mathieu@0: 	if (!time.empty()){
mathieu@0: #ifdef VERBOSE
mathieu@0: 		cout << "Time vector reset" << endl;
mathieu@0: #endif
mathieu@0: 		time.clear(); //clears the vector's elements
mathieu@0: 	}
mathieu@0: 	
mathieu@0: #ifdef VERBOSE
mathieu@0: 	cout << "End of get remaining features" << endl << endl;
mathieu@0: #endif
mathieu@0: 	return returnFeatures;
mathieu@0: }
mathieu@0: