annotate CalciumSignalAnalyser.cpp @ 4:f19dfc9cbd91 tip

Win32 build added

/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */

/*
 Calcium Signal Analyser Vamp Plugin
 
 Transient detection and characterisation for calcium signals.
 
 Centre for Digital Music, Queen Mary University of London.
 This file copyright 2010-2011 Mathieu Barthet and QMUL.
 
 Based on the QM Vamp note onset detector plugin by Christian 
 Landone, Chris Duxbury, and Juan Pablo Bello.
 
 This program 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 2 of the
 License, or (at your option) any later version.  See the file
 COPYING included with this distribution for more information.
  
 Version: 2
*/

#include "CalciumSignalAnalyser.h"

#include <dsp/onsets/PeakPicking.h>
#include <dsp/tempotracking/TempoTrack.h>

using std::string;
using std::vector;
using std::cerr;
using std::endl;

#define MAX_PEAKFREQ 10; //maximum peak frequency in Hz

//#define VERBOSE //remove comment to have output messages for data display or debugging information

CalciumSignalAnalyser::CalciumSignalAnalyser(float inputSampleRate) :
Vamp::Plugin(inputSampleRate),
m_inputSampleRate(inputSampleRate),
m_blockSize(0),
m_stepSize(0),
data(0),
time(0),
m_sensitivity(50),
m_delta(0.2),
m_mfwindur(10),
processcounter(0)
{
#ifdef VERBOSE
	cout << "Sample rate: " << inputSampleRate << endl;
#endif
}

CalciumSignalAnalyser::~CalciumSignalAnalyser()
{	
#ifdef VERBOSE
	cout << "Destructor" << endl;
#endif
	
	if (!data.empty()){
#ifdef VERBOSE
		cout << "Data vector reset" << endl;
#endif
		data.clear(); //clears the vector's elements
	}
	
	if (!time.empty()){
#ifdef VERBOSE
		cout << "Time vector reset" << endl;
#endif
		time.clear(); //clears the vector's elements
	}
}

string
CalciumSignalAnalyser::getIdentifier() const
{
    return "peakdetection";
}

string
CalciumSignalAnalyser::getName() const
{
    return "Calcium Signal Analyser";
}

string
CalciumSignalAnalyser::getDescription() const
{
    return "Automatic detection and characterisation of the signal's transients";
}

string
CalciumSignalAnalyser::getMaker() const
{
    return "Mathieu Barthet";
}

int
CalciumSignalAnalyser::getPluginVersion() const
{
    return 2;
}

string
CalciumSignalAnalyser::getCopyright() const
{
    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";
}

CalciumSignalAnalyser::ParameterList
CalciumSignalAnalyser::getParameterDescriptors() const
{
    ParameterList list;
	ParameterDescriptor desc;
	
    desc.identifier = "sensitivity";
    desc.name = "Peak-picking sensitivity";
    desc.description = "Sensitivity of the peak detection (criterion based on quadratic interpolation)";
    desc.minValue = 0;
    desc.maxValue = 100;
    desc.defaultValue = 50;
    desc.isQuantized = true;
    desc.quantizeStep = 1;
    desc.unit = "%";
    list.push_back(desc);
	
	desc.identifier = "deltathreshold";
    desc.name = "Peak-picking offset threshold (delta)";
    desc.description = "Offset threshold aiming at improving the peak detection by discarding the noise";
    desc.minValue = 0;
    desc.maxValue = 1;
    desc.defaultValue = 0.2;
    desc.isQuantized = true;
    desc.quantizeStep = 0.05;
    desc.unit = "";
    list.push_back(desc);
	
	desc.identifier = "mfwindowduration";
    desc.name = "Duration of the median filtering window";
    desc.description = "Sets the duration of the median filtering window";
    desc.minValue = 2;
    desc.maxValue = 100;
    desc.defaultValue = 10;
    desc.isQuantized = true;
    desc.quantizeStep = 0.1;
    desc.unit = "s";
    list.push_back(desc);

    return list;
}

float
CalciumSignalAnalyser::getParameter(std::string name) const
{	
	if (name == "sensitivity"){
		return m_sensitivity;
	}
	else if (name == "deltathreshold") {
		return m_delta;	
	}
	else if (name == "mfwindowduration") {
		return m_mfwindur;	
	}	
	
    return 0.0;
}

void
CalciumSignalAnalyser::setParameter(std::string name, float value)
{	
	if (name == "sensitivity") {
        if (m_sensitivity == value) return;
        m_sensitivity = value;
        
	}
	else if (name == "deltathreshold") {
		if (m_delta == value) return;
        m_delta = value;
	}
	else if (name == "mfwindowduration") {
		if (m_mfwindur == value) return;
        m_mfwindur = value;
	}
}


CalciumSignalAnalyser::ProgramList
CalciumSignalAnalyser::getPrograms() const
{
    ProgramList programs;

    return programs;
}

std::string
CalciumSignalAnalyser::getCurrentProgram() const
{	
	return "";
}

void
CalciumSignalAnalyser::selectProgram(std::string program)
{
}

bool
CalciumSignalAnalyser::initialise(size_t channels, size_t stepSize, size_t blockSize)
{
	
#ifdef VERBOSE
	cout << "Initialise" << endl;
#endif

	m_blockSize = 1;
	m_stepSize = 1;
	
#ifdef VERBOSE
	cout << m_blockSize << " " << m_stepSize << endl;
#endif
	
    if (channels < getMinChannelCount() ||
		channels > getMaxChannelCount()) {
        std::cerr << "CalciumSignalAnalyser::initialise: Unsupported channel count: "
		<< channels << std::endl;
        return false;
    }
	
    if (stepSize != getPreferredStepSize()) {
		
		std::cerr << "ERROR: CalciumSignalAnalyser::initialise: the step size has to be 1." << std::endl;
		return false;
	}
	
    if (blockSize != getPreferredBlockSize()) {
		
		std::cerr << "ERROR: CalciumSignalAnalyser::initialise: the block size has to be 1." << std::endl;
		return false;
	}
	
    return true;
}

void
CalciumSignalAnalyser::reset()
{	
#ifdef VERBOSE
	cout << "Reset" << endl;
#endif
	
	if (!data.empty()){
#ifdef VERBOSE
		cout << "Data vector reset" << endl;
#endif
		data.clear(); //clears the vector's elements
	}
	
	if (!time.empty()){
#ifdef VERBOSE
		cout << "Time vector reset" << endl;
#endif
		time.clear(); //clears the vector's elements
	}
}

size_t
CalciumSignalAnalyser::getPreferredStepSize() const
{	
	return getPreferredBlockSize();
}

size_t
CalciumSignalAnalyser::getPreferredBlockSize() const
{
	return 1;
}

size_t
CalciumSignalAnalyser::getMinChannelCount() const
{
    return 1;
}

size_t
CalciumSignalAnalyser::getMaxChannelCount() const
{
    return 1;
}


CalciumSignalAnalyser::OutputList
CalciumSignalAnalyser::getOutputDescriptors() const
{	
#ifdef VERBOSE
	cout << "get output desc" << endl;
#endif
	
    OutputList list;
  	
    OutputDescriptor pt;
    pt.identifier = "peaktimes";
    pt.name = "Peaks times";
    pt.description = "Time positions of the peaks (in seconds)";
    pt.unit = "";
    pt.hasFixedBinCount = true;
    pt.binCount = 0; //0 for time only feature
    pt.sampleType = OutputDescriptor::VariableSampleRate;
	pt.hasDuration = false;
	pt.hasKnownExtents = false;
	pt.isQuantized = false;
    pt.sampleRate = 0; //no sensible resolution, the time of the feature are given in the timestamp of the feature
	list.push_back(pt);
	
	OutputDescriptor po;
    po.identifier = "peakonsets";
    po.name = "Peaks onsets";
    po.description = "Time positions of the peak onsets (in seconds)";
    po.unit = "";
    po.hasFixedBinCount = true;
    po.binCount = 0; //0 for time only feature
    po.sampleType = OutputDescriptor::VariableSampleRate;
	po.hasDuration = false;
	po.hasKnownExtents = false;
	po.isQuantized = false;
    po.sampleRate = 0; //no sensible resolution, the time of the feature are given in the timestamp of the feature
	list.push_back(po);
	
	OutputDescriptor sdf;
    sdf.identifier = "smoothed_df";
    sdf.name = "Smoothed Detection Function";
    sdf.description = "Smoothed probability function used for peak-picking";
    sdf.unit = "";
    sdf.hasFixedBinCount = true;
    sdf.binCount = 1;
    sdf.hasKnownExtents = false;
    sdf.isQuantized = false;
    sdf.sampleType = OutputDescriptor::FixedSampleRate;
	sdf.sampleRate = m_inputSampleRate; //the smoothed detection function has the same sample rate as the input signal
    list.push_back(sdf);
		
	OutputDescriptor pf;
    pf.identifier = "peak_frequency";
    pf.name = "Peak frequency";
    pf.description = "Frequency of the peaks (in Hz)";
    //pf.unit = "Hz"; //if unit is Hz then the feature is treated as a pitch which is not relevant in this case
    pf.unit = "";
	pf.hasFixedBinCount = true;
    pf.binCount = 1; //1 corresponds to time and value
    pf.sampleType = OutputDescriptor::VariableSampleRate;
	pf.sampleRate = 0; //no sensible resolution, the time of the feature is given in the timestamp of the feature
	pf.hasKnownExtents = true;
	pf.minValue = 0;
	pf.maxValue = MAX_PEAKFREQ;	
	pf.isQuantized = false;
	pf.hasDuration = true; //duration of the signal
	list.push_back(pf);
    
	OutputDescriptor mipi;
    mipi.identifier = "mean_interpeakinterval";
    mipi.name = "Mean Inter Peak Interval";
    mipi.description = "Mean of the time intervals between peaks (in seconds)";
    mipi.unit = "s";
	mipi.hasFixedBinCount = true;
    mipi.binCount = 1; //1 corresponds to time and value
    mipi.sampleType = OutputDescriptor::VariableSampleRate;
	mipi.sampleRate = 0; //no sensible resolution, the time of the feature is given in the timestamp of the feature
	mipi.hasKnownExtents = false;
	mipi.isQuantized = false;
	mipi.hasDuration = true; //duration of the signal
	list.push_back(mipi);
	
    return list;
}

CalciumSignalAnalyser::FeatureSet
CalciumSignalAnalyser::process(const float *const *inputBuffers,
                       Vamp::RealTime timestamp)
{	
	FeatureSet fs;
	
	if (m_blockSize == 0 || m_stepSize == 0) {
		cerr << "ERROR: CalciumSignalAnalyser::process: "
		<< "CalciumSignalAnalyser has not been initialised."
		<< endl;
		return fs;
    }
	
	data.push_back(inputBuffers[0][0]); //assumes that only one sample is read at a time
	time.push_back(timestamp);
	
	processcounter ++;
	
	fs = FeatureSet(); //empty feature set
	
	return fs;	
}

CalciumSignalAnalyser::FeatureSet
CalciumSignalAnalyser::getRemainingFeatures()
{
   	FeatureSet returnFeatures;
	returnFeatures = FeatureSet(); 
	
#ifdef VERBOSE
	cout << "Get remaining features" << endl << endl;
#endif
	
	/***************************************
	 Peak times
	 ****************************************/
    
	double aCoeffs[] = { 1.0000, -0.5949, 0.2348 };
    double bCoeffs[] = { 0.1600,  0.3200, 0.1600 };
	
    //The signal is directly treated as a detection function for the peak picking stage
	
    PPickParams ppParams;
    ppParams.length = data.size(); //length of the signal
	
#ifdef VERBOSE
	cout << "Length: " << ppParams.length << endl;
#endif
	
	// tau and cutoff appear to be unused in PeakPicking, but I've
    // inserted some moderately plausible values rather than leave
    // them unset.  The QuadThresh values come from trial and error.
    // The rest of these are copied from ttParams in the BeatTracker
    // code: I don't claim to know whether they're good or not --cc
        
	ppParams.tau = m_stepSize / m_inputSampleRate; //not used for peak picking 
	
	ppParams.alpha = 9;
    ppParams.cutoff = m_inputSampleRate/4; //not used for peak picking 
    ppParams.LPOrd = 2;
    ppParams.LPACoeffs = aCoeffs;
    ppParams.LPBCoeffs = bCoeffs;
    
	int NWin;
	float samples = m_inputSampleRate*m_mfwindur;
	if (samples < 4){
		NWin = 4; //minimal possible length of the window in samples
#ifdef VERBOSE
		cout << "WARNING: CalciumSignalAnalyser::getRemainingFeatures:"
		<< "The duration of the median filtering window is too small considering the"
		<< "sample frequency of the signal. Duration of "<< NWin/m_inputSampleRate
		<< " s used instead."
		<< endl;
#endif
	}else{
		NWin = floor(m_inputSampleRate*m_mfwindur);
		if (NWin%2 != 0) NWin--;
	}	
		
	ppParams.WinT.post = NWin/2;
    ppParams.WinT.pre = NWin/2 - 1;
    		
	ppParams.QuadThresh.a = (100 - m_sensitivity) / 1000.0;
    ppParams.QuadThresh.b = 0;
    ppParams.QuadThresh.c = (100 - m_sensitivity) / 1500.0;
	ppParams.delta = m_delta; //delta threshold used as an offset when computing the smoothed detection function
			
    PeakPicking peakPicker(ppParams);
	
	double *ppSrc = new double[ppParams.length];
	if (ppSrc==NULL)	{
		cerr << "ERROR: CalciumSignalAnalyser::getRemainingFeatures:"
		<< "Unable to allocate memory (ppSrc)."
		<< endl;
		return returnFeatures;
	}
		
    for (unsigned int i = 0; i < ppParams.length; ++i) {
        
		ppSrc[i] = (double) data[i];
	}	
	
    vector<int> peaks;
	vector<int> peak_frame;
	vector<int> onsets_frame;
	
	peakPicker.process(ppSrc, ppParams.length, peaks);
	
	for (size_t i = 0; i < peaks.size(); ++i) {
		
		//peak times
        size_t index = peaks[i];
		size_t peakframe = peaks[i] * m_stepSize;
		peak_frame.push_back(peakframe);
		
		Feature featurepeaks;
		featurepeaks.hasTimestamp = true;	
		featurepeaks.timestamp = Vamp::RealTime::frame2RealTime(peakframe, lrintf(m_inputSampleRate));	
		
		returnFeatures[0].push_back(featurepeaks); // peak times are output 0
		
		//algorithm to detect the peak onsets
		 
		double prevDiff = 0.0;
		while (index > 1) {
			double diff = ppSrc[index] - ppSrc[index-1];
			if (diff < prevDiff * 0.9) break;
			prevDiff = diff;
			--index;
		}	
			
		size_t onsetframe = index * m_stepSize;		
		onsets_frame.push_back(onsetframe);
		
		Feature featureonsets;
		featureonsets.hasTimestamp = true;	
		featureonsets.timestamp = Vamp::RealTime::frame2RealTime(onsetframe, lrintf(m_inputSampleRate));	
		
#ifdef VERBOSE
		cout << featureonsets.timestamp << endl;
#endif
		returnFeatures[1].push_back(featureonsets); // peak onsets are output 1
    }
	
    for (unsigned int i = 0; i < ppParams.length; ++i) {
        
        Feature feature;
		feature.hasTimestamp = true;
		size_t frame = i * m_stepSize;
				
		feature.timestamp = Vamp::RealTime::frame2RealTime(frame, lrintf(m_inputSampleRate));
		
        feature.values.push_back(ppSrc[i]);
     	
		returnFeatures[2].push_back(feature); // smoothed df is output 2
    }
	
	/***************************************
	 Peak frequency
	 ****************************************/
	
	double peakfreq, duration;
	int npeaks; //total number of peaks
	
	//last element of the time vector is the duration of the signal - converts the RealTime to a float
	duration = time.back()/Vamp::RealTime(1,0); //division by 1 s to obtain a double
		
	npeaks = peaks.size();
	peakfreq = npeaks/duration;
	
#ifdef VERBOSE
	cout << npeaks << endl;
	cout << duration << endl;	
	cout << peakfreq << endl;
#endif
	
	Feature pf;
	pf.hasTimestamp = true;
	pf.timestamp = time[0];
	pf.hasDuration = true;			
	pf.duration = time.back()-time[0];
	
#ifdef VERBOSE
	cout << "processcounter" << processcounter << endl;
	
	cout << "time.back " << time.back().toText() << endl;
	cout << "time.size " << time.size() << endl;
	cout << "data.size " << data.size() << endl;
	cout << "time[0, 1 ...] " << time[0].toText() << " " << time[1].toText() << " " << time[2].toText() << " " << time[time.size()].toText() << endl;
	
	cout << "time[999] " << time[999].toText() << endl;
	cout << "data[999] " << data[999] << endl;
	cout << "time[1000] " << time[1000].toText() << endl;
	cout << "data[1000] " << data[1000]<< endl;
#endif
	
	pf.values.push_back(peakfreq);
		
	char pflabel[256];
	sprintf(pflabel,"%f Hz",peakfreq);
	pf.label = pflabel;
	
#ifdef VERBOSE
	cout << "label: " << pflabel << endl;
#endif
	
	returnFeatures[3].push_back(pf); //peak frequency is output 3
	
	/***************************************
	 Mean Inter Peak Interval
	 ****************************************/
	
	double o1,o2;
	
	double sumipi_s = 0;
	double meanipi_s;
	
	for (int j = 0; j < npeaks-1; ++j) {
		o2 = Vamp::RealTime::frame2RealTime(peak_frame[j+1], lrintf(m_inputSampleRate))/Vamp::RealTime(1,0);
		o1 = Vamp::RealTime::frame2RealTime(peak_frame[j], lrintf(m_inputSampleRate))/Vamp::RealTime(1,0);
#ifdef VERBOSE
		cout << "o1: " << o1 << endl;
		cout << "o2: " << o2 << endl;
#endif
		sumipi_s = sumipi_s + o2 - o1;	
	}
	
	meanipi_s = sumipi_s/(npeaks-1);
	
#ifdef VERBOSE
	cout << "Sum IPI: " << sumipi_s << endl;
	cout << "Mean IPI: " << meanipi_s << endl;
#endif
	
	Feature mipi;
	mipi.hasTimestamp = true;
	mipi.timestamp = time[0];
	mipi.hasDuration = true;			
	mipi.duration = time.back()-time[0];
	
	mipi.values.push_back(meanipi_s); //mean inter peak interval in secs
	
	returnFeatures[4].push_back(mipi); //mean inter peak interval is output 4
	
	//clear data
	
	delete [] ppSrc;
	
	if (!data.empty()){
#ifdef VERBOSE
		cout << "Data vector reset" << endl;
#endif
		data.clear(); //clears the vector's elements
	}
	
	if (!time.empty()){
#ifdef VERBOSE
		cout << "Time vector reset" << endl;
#endif
		time.clear(); //clears the vector's elements
	}
	
#ifdef VERBOSE
	cout << "End of get remaining features" << endl << endl;
#endif
	return returnFeatures;
}