Skip to content

GitLab

Commit d438192b authored by nfoy's avatar nfoy
Browse files

Merge branch 'wip-jlindgre-merge-fastica-mods' into integration-1.2

parents 9dd3b1c8 5855ad3f
Expand all Hide whitespace changes
Inline Side-by-side
contrib/plugins/processing/signal-processing/doc/Doc_BoxAlgorithm_IndependentComponentAnalysisFastICA.dox-part
View file @ d438192b
/**
* \page BoxAlgorithm_IndependentComponentAnalysisFastICA Independent Component AnalysisFast (FastICA) -
* \page BoxAlgorithm_IndependentComponentAnalysisFastICA Independent Component Analysis (FastICA)
__________________________________________________________________
Detailed description
__________________________________________________________________
* |OVP_DocBegin_BoxAlgorithm_IndependentComponentAnalysisFastICA_Description|
* This plugin is used to decompose into independent components the
* input signal. This plugin is based on the Fast ICA algorithm.
* This box attempts to find a decomposition of the signal to its
* independent components. The approach is based on the FastICA algorithm.
* |OVP_DocEnd_BoxAlgorithm_IndependentComponentAnalysisFastICA_Description|
__________________________________________________________________
......@@ -29,6 +29,7 @@ __________________________________________________________________
* |OVP_DocEnd_BoxAlgorithm_IndependentComponentAnalysisFastICA_Outputs|
*
* |OVP_DocBegin_BoxAlgorithm_IndependentComponentAnalysisFastICA_Output1|
* The decomposed signal.
* |OVP_DocEnd_BoxAlgorithm_IndependentComponentAnalysisFastICA_Output1|
__________________________________________________________________
......@@ -36,16 +37,68 @@ Settings description
__________________________________________________________________
* |OVP_DocBegin_BoxAlgorithm_IndependentComponentAnalysisFastICA_Settings|
*
* |OVP_DocBegin_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting1|
* Number of independent components to extract (equals PCA dimension reduction)
* |OVP_DocEnd_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting1|
*
* |OVP_DocBegin_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting2|
* How many seconds of sample to collect to estimate the ICA model?
* |OVP_DocEnd_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting2|
*
* |OVP_DocBegin_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting3|
* Decomposition type. Deflation is an approach where each component is
* estimated separately in turns. Symmetric estimation optimizes all
* components at once.
* |OVP_DocEnd_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting3|
*
* |OVP_DocBegin_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting4|
* Maximum number of iterations
* |OVP_DocEnd_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting4|
*
* |OVP_DocBegin_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting5|
* Enable fine tuning?
* |OVP_DocEnd_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting5|
*
* |OVP_DocBegin_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting6|
* Maximum number of iterations for the fine tuning
* |OVP_DocEnd_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting6|
*
* |OVP_DocBegin_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting7|
* Used nonlinearity type
* |OVP_DocEnd_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting7|
*
* |OVP_DocBegin_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting8|
* Mu parameter
* |OVP_DocEnd_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting8|
*
* |OVP_DocBegin_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting9|
* Epsilon parameter
* |OVP_DocEnd_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting9|
*
* |OVP_DocBegin_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting10|
* Filename to save the estimated decomposition matrix W to
* |OVP_DocEnd_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting10|
*
* |OVP_DocBegin_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting11|
* Should the matrix W be saved to a file?
* |OVP_DocEnd_BoxAlgorithm_IndependentComponentAnalysisFastICA_Setting11|
*
* |OVP_DocEnd_BoxAlgorithm_IndependentComponentAnalysisFastICA_Settings|
__________________________________________________________________
Examples description
__________________________________________________________________
* |OVP_DocBegin_BoxAlgorithm_IndependentComponentAnalysisFastICA_Examples|
* Let's consider our input signal.
* To decompose the input signal and to enlight for example the blinks,
* apply the FastICA. One of the new component will contain blinks.
* One use-case of ICA is to attempt to separate the signal of interest from
* nuisance artifacts. For example, supposing that ICA makes a meaningful decomposition
* of your EEG signal, you will see artifacts such as those from eyeblinks more
* clearly segregated to specific output channels instead of contaminating
* all of the channels.
* |OVP_DocEnd_BoxAlgorithm_IndependentComponentAnalysisFastICA_Examples|
__________________________________________________________________
......@@ -53,6 +106,20 @@ Miscellaneous description
__________________________________________________________________
* |OVP_DocBegin_BoxAlgorithm_IndependentComponentAnalysisFastICA_Miscellaneous|
* This plugin applies the Fast ICA algorithm to the input signal.
* This plugin applies the FastICA algorithm to the input signal. The box can store the
* estimated decomposition matrix W to a file. This file can then be used later
* in the spatial filter box to apply the decomposition on fresh data.
*
* The box also outputs the decomposed signal, but the decomposition is active only
* after the model has been estimated (after the specified number of samples have been collected).
* If you wish to decompose the whole data, then you can first train the ICA model, save the matrix,
* and then separately apply it to the original data with the spatial filter.
*
* The FastICA algorithm is described in
*
* A. Hyvrinen. "Fast and Robust Fixed-Point Algorithms for Independent Component Analysis", IEEE Transactions on Neural Networks 10(3):626-634, 1999.
*
* The implementation used by the box is from the ITPP toolkit.
*
* |OVP_DocEnd_BoxAlgorithm_IndependentComponentAnalysisFastICA_Miscellaneous|
*/
contrib/plugins/processing/signal-processing/src/box-algorithms/ovpCFastICA.cpp
View file @ d438192b
......@@ -12,6 +12,7 @@ using namespace itpp;
using namespace OpenViBE;
using namespace OpenViBE::Plugins;
using namespace OpenViBE::Kernel;
using namespace OpenViBEPlugins;
using namespace OpenViBEPlugins::SignalProcessing;
......@@ -23,27 +24,90 @@ void CFastICA::computeICA(void)
const uint32 l_ui32ChannelCount = m_oDecoder.getOutputMatrix()->getDimensionSize(0);
const uint32 l_ui32SampleCount = m_oDecoder.getOutputMatrix()->getDimensionSize(1);
const float64 *l_pInputBuffer = m_oDecoder.getOutputMatrix()->getBuffer();
float64 *l_pOutputBuffer = m_oEncoder.getInputMatrix()->getBuffer();
mat sources(l_ui32ChannelCount, l_ui32SampleCount);
mat ICs(l_ui32ChannelCount, l_ui32SampleCount);
const uint32 l_ui32NumOfICs = m_oEncoder.getInputMatrix()->getDimensionSize(0);
mat sources(l_ui32ChannelCount, l_ui32SampleCount); // current block (for decomposing)
mat Buffer_sources(l_ui32ChannelCount, m_ui32Buff_Size); // accumulated blocks (for training)
mat ICs(l_ui32NumOfICs, l_ui32SampleCount);
//mat Mix_mat(l_ui32ChannelCount, l_ui32ChannelCount);
mat Sep_mat(l_ui32NumOfICs, l_ui32ChannelCount);
//mat Dewhite(l_ui32ChannelCount, l_ui32ChannelCount);
// Append the data to a FIFO buffer
for (uint32 i=0; i < l_ui32ChannelCount; i++)
{
for(uint32 j=0 ; j<l_ui32SampleCount ; j++)
for(uint32 j=0 ; j<m_ui32Buff_Size ; j++)
{
sources((int)i, (int)j) = (double)l_pInputBuffer[i*l_ui32SampleCount+j];
if(j<m_ui32Buff_Size-l_ui32SampleCount)
{
m_pFifoBuffer[i*m_ui32Buff_Size+j+l_ui32SampleCount]=m_pFifoBuffer[i*m_ui32Buff_Size+j]; // memory shift
if(j<l_ui32SampleCount)
{
m_pFifoBuffer[i*m_ui32Buff_Size+j] = (double)l_pInputBuffer[i*l_ui32SampleCount+l_ui32SampleCount-1-j];
sources((int)i, (int)j) = (double)l_pInputBuffer[i*l_ui32SampleCount+j];
}
}
Buffer_sources((int)i, (int)(m_ui32Buff_Size-1-j)) = m_pFifoBuffer[i*m_ui32Buff_Size+j];
}
}
Fast_ICA fastica(sources);
fastica.set_nrof_independent_components(sources.rows());
fastica.set_non_linearity(FICA_NONLIN_TANH);
fastica.set_approach(FICA_APPROACH_DEFL);
fastica.separate();
ICs = fastica.get_independent_components();
m_ui32Samp_Nb += l_ui32SampleCount;
if((m_ui32Samp_Nb >= m_ui32Buff_Size)&&(m_bTrained==false))
{
this->getLogManager() << LogLevel_Trace << "Instanciating the Fast_ICA object with " << m_ui32Samp_Nb << " samples.\n";
Fast_ICA fastica(Buffer_sources);
this->getLogManager() << LogLevel_Trace << "Setting the number of ICs to extract to " << l_ui32NumOfICs << " and configuring FastICA...\n";
fastica.set_nrof_independent_components(l_ui32NumOfICs);
fastica.set_approach(m_ui32Type);
fastica.set_non_linearity(m_ui32Non_Lin);
fastica.set_max_num_iterations(m_ui32NbRep_max);
fastica.set_fine_tune(m_bSetFineTune);
fastica.set_max_fine_tune(m_ui32NbTune_max);
fastica.set_mu(m_ui64Set_Mu);
fastica.set_epsilon(m_ui64Epsilon);
//if(m_ui32Samp_Nb>l_ui32SampleCount) fastica.set_init_guess((Dewhite * Dewhite.T()) * Sep_mat.T());
this->getLogManager() << LogLevel_Trace << "Explicit launch of the Fast_ICA algorithm. Can occasionnally take time.\n";
fastica.separate();
this->getLogManager() << LogLevel_Trace << "Retrieving ICs from fastica .\n";
//ICs = fastica.get_independent_components();
//this->getLogManager() << LogLevel_Trace << "Retrieving mixing matrix from fastica .\n";
//Mix_mat = fastica.get_mixing_matrix();
this->getLogManager() << LogLevel_Trace << "Retrieving separating matrix from fastica .\n";
Sep_mat = fastica.get_separating_matrix();
//Dewhite = fastica.get_dewhitening_matrix();
m_bTrained = true;
float64 *l_pDemixer = m_oDemixer.getBuffer();
for(uint32 i=0; i < l_ui32NumOfICs; i++)
{
for(uint32 j=0; j < l_ui32ChannelCount; j++)
{
l_pDemixer[i*l_ui32ChannelCount+j] = Sep_mat((int)i,(int)j);
}
}
}
else
{
// Use the previously stored matrix
const float64* l_pDemixer = m_oDemixer.getBuffer();
for(uint32 i=0;i<l_ui32NumOfICs;i++)
{
for(uint32 j=0;j<l_ui32ChannelCount;j++)
{
Sep_mat((int)i, (int)j) = l_pDemixer[i*l_ui32ChannelCount+j];
}
}
}
for (uint32 i=0; i < l_ui32ChannelCount; i++)
// Effective demixing (ICA after m_ui32Duration sec)
ICs = Sep_mat * sources;
float64 *l_pOutputBuffer = m_oEncoder.getInputMatrix()->getBuffer();
//this->getLogManager() << LogLevel_Trace << "Filling output buffer with ICs .\n";
for (uint32 i=0; i < l_ui32NumOfICs; i++)
{
for(uint32 j=0 ; j < l_ui32SampleCount ; j++)
{
......@@ -66,6 +130,28 @@ boolean CFastICA::initialize()
m_oDecoder.initialize(*this, 0);
m_oEncoder.initialize(*this, 0);
m_ui32Nb_ICs = FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 0);
m_ui32Duration = FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 1);
uint32 l_ui32Type = FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 2);
m_ui32NbRep_max = FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 3);
m_bSetFineTune = FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 4);
m_ui32NbTune_max = FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 5);
m_ui32Non_Lin = FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 6);
m_ui64Set_Mu = FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 7);
m_ui64Epsilon = FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 8);
m_sSpatialFilterFilename = FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 9);
m_bSaveAsFile = FSettingValueAutoCast(*this->getBoxAlgorithmContext(), 10);
m_ui32Type = (l_ui32Type==0 ? FICA_APPROACH_DEFL : FICA_APPROACH_SYMM);
m_pFifoBuffer = NULL;
if(m_bSaveAsFile && m_sSpatialFilterFilename==CString(""))
{
this->getLogManager() << "If save is enabled, filename must be provided\n";
return false;
}
return true;
}
......@@ -73,6 +159,11 @@ boolean CFastICA::uninitialize()
{
m_oEncoder.uninitialize();
m_oDecoder.uninitialize();
if(m_pFifoBuffer)
{
delete[] m_pFifoBuffer;
m_pFifoBuffer = NULL;
}
return true;
}
......@@ -95,30 +186,83 @@ boolean CFastICA::process()
if(m_oDecoder.isHeaderReceived())
{
// Set the output (encoder) matrix prorperties from the input (decoder)
if(m_oDecoder.getOutputMatrix()->getDimensionCount()!=2)
{
this->getLogManager() << LogLevel_Error << "Needs a 2 dimensional (rows x cols) matrix as input\n";
return false;
}
m_ui32Buff_Size = static_cast<uint32>(m_oDecoder.getOutputSamplingRate())*m_ui32Duration;
const uint32 l_ui32ChannelCount = m_oDecoder.getOutputMatrix()->getDimensionSize(0);
const uint32 l_ui32SampleCount = m_oDecoder.getOutputMatrix()->getDimensionSize(1);
if(m_pFifoBuffer)
{
delete[] m_pFifoBuffer;
}
m_pFifoBuffer = new OpenViBE::float64[l_ui32ChannelCount*m_ui32Buff_Size];
this->getLogManager() << LogLevel_Trace << "FIFO buffer initialized with " << l_ui32ChannelCount*m_ui32Buff_Size << ".\n";
if(m_ui32Nb_ICs > l_ui32ChannelCount)
{
this->getLogManager() << LogLevel_Warning << "Trying to estimate more components than channels, truncating\n";
m_ui32Nb_ICs = l_ui32ChannelCount;
}
for(uint32 j=0 ; j<l_ui32ChannelCount*m_ui32Buff_Size ; j++)
{
m_pFifoBuffer[j]=0.0;
}
m_ui32Samp_Nb = 0;
m_bTrained = false;
IMatrix* l_pEncoderMatrix = m_oEncoder.getInputMatrix();
l_pEncoderMatrix->setDimensionCount(2);
l_pEncoderMatrix->setDimensionSize(0, m_ui32Nb_ICs);
l_pEncoderMatrix->setDimensionSize(1, l_ui32SampleCount);
OpenViBEToolkit::Tools::Matrix::copyDescription(*l_pEncoderMatrix, *m_oDecoder.getOutputMatrix());
m_oEncoder.getInputSamplingRate() = m_oDecoder.getOutputSamplingRate();
for(uint32 i=0 ; i<l_pEncoderMatrix->getDimensionSize(0) ; i++)
for(uint32 c=0 ; c < m_ui32Nb_ICs ; c++)
{
char l_sBuffer[64];
sprintf(l_sBuffer, "IC %d", i+1);
l_pEncoderMatrix->setDimensionLabel(0,i, l_sBuffer);
sprintf(l_sBuffer, "IC %d", c+1);
l_pEncoderMatrix->setDimensionLabel(0,c, l_sBuffer);
}
m_oEncoder.encodeHeader();
m_oDemixer.setDimensionCount(2);
m_oDemixer.setDimensionSize(0,m_ui32Nb_ICs);
m_oDemixer.setDimensionSize(1,l_ui32ChannelCount);
// Set the demixer to (partial) identity matrix to start with
OpenViBEToolkit::Tools::Matrix::clearContent(m_oDemixer);
float64* l_pDemixer = m_oDemixer.getBuffer();
for(uint32 c=0;c<m_ui32Nb_ICs;c++)
{
l_pDemixer[c*l_ui32ChannelCount+c] = 1.0;
}
getBoxAlgorithmContext()->getDynamicBoxContext()->markOutputAsReadyToSend(0, 0, 0);
}
if(m_oDecoder.isBufferReceived())
{
const uint64 l_ui64LastChunkStartTime = l_pDynamicBoxContext->getInputChunkStartTime(0,i);
const uint64 l_ui64LastChunkEndTime = l_pDynamicBoxContext->getInputChunkEndTime(0,i);
const uint64 l_ui64LastChunkEndTime = l_pDynamicBoxContext->getInputChunkEndTime(0,i);
const uint32 l_ui32ChannelCount = m_oDecoder.getOutputMatrix()->getDimensionSize(0);
computeICA();
if((m_bSaveAsFile) && (m_bTrained) && (m_bFileSaved==false))
{
if(!OpenViBEToolkit::Tools::Matrix::saveToTextFile(m_oDemixer, m_sSpatialFilterFilename))
this->getLogManager() << LogLevel_Warning << "Unable to save to [" << m_sSpatialFilterFilename << "\n";
m_bFileSaved=true;
}
m_oEncoder.encodeBuffer();
getBoxAlgorithmContext()->getDynamicBoxContext()->markOutputAsReadyToSend(0, l_ui64LastChunkStartTime, l_ui64LastChunkEndTime);
......
contrib/plugins/processing/signal-processing/src/box-algorithms/ovpCFastICA.h
View file @ d438192b
......@@ -7,13 +7,18 @@
#if defined TARGET_HAS_ThirdPartyITPP
#include "../ovp_defines.h"
#include <openvibe/ov_all.h>
#include <toolkit/ovtk_all.h>
#include <vector>
#include <map>
#include <string>
#define FICA_NONLIN_POW3 10 // Use x^3 non-linearity.
#define FICA_NONLIN_TANH 20 // Use tanh(x) non-linearity.
#define FICA_NONLIN_GAUSS 30 // Use Gaussian non-linearity.
#define FICA_NONLIN_SKEW 40 // Use skew non-linearity.
// TODO create a member function to get rid of this
#ifndef CString2Boolean
#define CString2Boolean(string) (strcmp(string,"true"))?0:1
......@@ -43,14 +48,34 @@ namespace OpenViBEPlugins
_IsDerivedFromClass_Final_(OpenViBE::Plugins::IBoxAlgorithm, OVP_ClassId_FastICA)
public:
protected:
virtual void computeICA(void);
public:
OpenViBEToolkit::TSignalDecoder<CFastICA> m_oDecoder;
OpenViBEToolkit::TSignalEncoder<CFastICA> m_oEncoder;
protected:
OpenViBE::float64* m_pFifoBuffer;
OpenViBE::CMatrix m_oDemixer; // The estimated matrix W
bool m_bTrained;
bool m_bFileSaved;
OpenViBE::uint32 m_ui32Buff_Size;
OpenViBE::uint32 m_ui32Samp_Nb;
OpenViBE::uint32 m_ui32Nb_ICs;
OpenViBE::uint32 m_ui32Duration;
OpenViBE::uint32 m_ui32NbRep_max;
OpenViBE::uint32 m_ui32NbTune_max;
OpenViBE::CString m_sSpatialFilterFilename;
OpenViBE::boolean m_bSaveAsFile;
OpenViBE::boolean m_bSetFineTune;
OpenViBE::float64 m_ui64Set_Mu;
OpenViBE::float64 m_ui64Epsilon;
OpenViBE::uint32 m_ui32Non_Lin;
OpenViBE::uint32 m_ui32Type;
};
......@@ -59,13 +84,13 @@ namespace OpenViBEPlugins
public:
virtual void release(void) { }
virtual OpenViBE::CString getName(void) const { return OpenViBE::CString("Independent component analysis (FastICA)"); }
virtual OpenViBE::CString getAuthorName(void) const { return OpenViBE::CString("Guillaume Gibert"); }
virtual OpenViBE::CString getAuthorCompanyName(void) const { return OpenViBE::CString("INSERM"); }
virtual OpenViBE::CString getName(void) const { return OpenViBE::CString("Independent Component Analysis (FastICA)"); }
virtual OpenViBE::CString getAuthorName(void) const { return OpenViBE::CString("Guillaume Gibert / Jeff B."); }
virtual OpenViBE::CString getAuthorCompanyName(void) const { return OpenViBE::CString("INSERM / Independent"); }
virtual OpenViBE::CString getShortDescription(void) const { return OpenViBE::CString("Computes fast independent component analysis"); }
virtual OpenViBE::CString getDetailedDescription(void) const { return OpenViBE::CString(""); }
virtual OpenViBE::CString getCategory(void) const { return OpenViBE::CString("Signal processing/Independent component analysis"); }
virtual OpenViBE::CString getVersion(void) const { return OpenViBE::CString("0.1"); }
virtual OpenViBE::CString getVersion(void) const { return OpenViBE::CString("0.2"); }
virtual OpenViBE::CIdentifier getCreatedClass(void) const { return OVP_ClassId_FastICA; }
virtual OpenViBE::Plugins::IPluginObject* create(void) { return new OpenViBEPlugins::SignalProcessing::CFastICA(); }
......@@ -74,6 +99,19 @@ namespace OpenViBEPlugins
{
rPrototype.addInput ("Input signal", OV_TypeId_Signal);
rPrototype.addOutput("Output signal", OV_TypeId_Signal);
rPrototype.addSetting("Number of independent components to extract", OV_TypeId_Integer, "14");
rPrototype.addSetting("Sample size (seconds) for estimation", OV_TypeId_Integer, "120");
rPrototype.addSetting("Decomposition type (0==deflate, 1==symmetric)", OV_TypeId_Integer, "1");
rPrototype.addSetting("Max number of reps for the ICA convergence", OV_TypeId_Integer, "100000");
rPrototype.addSetting("Fine tuning", OV_TypeId_Boolean, "true");
rPrototype.addSetting("Max number of reps for the fine tuning", OV_TypeId_Integer, "100");
rPrototype.addSetting("Non linearity (10: POW3, 20: TANH, 30: GAUSS)", OV_TypeId_Integer, "20");
rPrototype.addSetting("Internal Mu parameter for FastICA", OV_TypeId_Float, "1.0");
rPrototype.addSetting("Internal Epsilon parameter for FastICA", OV_TypeId_Float, "0.0001");
rPrototype.addSetting("Spatial filter filename", OV_TypeId_Filename, "");
rPrototype.addSetting("Save the spatial filter/demixing matrix", OV_TypeId_Boolean, "true");
rPrototype.addFlag (OpenViBE::Kernel::BoxFlag_IsUnstable);
return true;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment