ovpCHilbertTransform.cpp 5.06 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147
#if defined(TARGET_HAS_ThirdPartyEIGEN)

#include "ovpCHilbertTransform.h"
#include <complex>
#include <Eigen/Dense>
#include <unsupported/Eigen/FFT>
#include <iostream>

using namespace OpenViBE;
using namespace OpenViBE::Kernel;
using namespace OpenViBE::Plugins;

using namespace OpenViBEPlugins;
using namespace OpenViBEPlugins::SignalProcessing;

using namespace Eigen;


boolean CAlgorithmHilbertTransform::initialize(void)
{
	ip_pMatrix.initialize(this->getInputParameter(OVP_Algorithm_HilbertTransform_InputParameterId_Matrix));
	op_pHilbertMatrix.initialize(this->getOutputParameter(OVP_Algorithm_HilbertTransform_OutputParameterId_HilbertMatrix));
	op_pEnvelopeMatrix.initialize(this->getOutputParameter(OVP_Algorithm_HilbertTransform_OutputParameterId_EnvelopeMatrix));
	op_pPhaseMatrix.initialize(this->getOutputParameter(OVP_Algorithm_HilbertTransform_OutputParameterId_PhaseMatrix));


	return true;
}

boolean CAlgorithmHilbertTransform::uninitialize(void)
{
	op_pHilbertMatrix.uninitialize();
	op_pEnvelopeMatrix.uninitialize();
	op_pPhaseMatrix.uninitialize();
	ip_pMatrix.uninitialize();
	return true;
}

boolean CAlgorithmHilbertTransform::process(void)
{

	uint32 l_ui32ChannelCount = ip_pMatrix->getDimensionSize(0);
	uint32 l_ui32SamplesPerChannel = ip_pMatrix->getDimensionSize(1);


	IMatrix* l_pInputMatrix = ip_pMatrix;
	IMatrix* l_pOutputHilbertMatrix = op_pHilbertMatrix;
	IMatrix* l_pOutputEnvelopeMatrix = op_pEnvelopeMatrix;
	IMatrix* l_pOutputPhaseMatrix = op_pPhaseMatrix;

	if(this->isInputTriggerActive(OVP_Algorithm_HilbertTransform_InputTriggerId_Initialize)) //Check if the input is correct
	{
		if( l_pInputMatrix->getDimensionCount() != 2)
		{
			this->getLogManager() << LogLevel_Error << "The input matrix must have 2 dimensions, here the dimension is "<<l_pInputMatrix->getDimensionCount()<<"\n";
			return false;
		}

		//Setting size of outputs

		l_pOutputHilbertMatrix->setDimensionCount(2);
		l_pOutputHilbertMatrix->setDimensionSize(0,l_ui32ChannelCount);
		l_pOutputHilbertMatrix->setDimensionSize(1,l_ui32SamplesPerChannel);

		l_pOutputEnvelopeMatrix->setDimensionCount(2);
		l_pOutputEnvelopeMatrix->setDimensionSize(0,l_ui32ChannelCount);
		l_pOutputEnvelopeMatrix->setDimensionSize(1,l_ui32SamplesPerChannel);

		l_pOutputPhaseMatrix->setDimensionCount(2);
		l_pOutputPhaseMatrix->setDimensionSize(0,l_ui32ChannelCount);
		l_pOutputPhaseMatrix->setDimensionSize(1,l_ui32SamplesPerChannel);

		for(uint32 i=0; i<l_ui32ChannelCount;i++)
		{
			l_pOutputHilbertMatrix->setDimensionLabel(0,i,l_pInputMatrix->getDimensionLabel(0,i));
			l_pOutputEnvelopeMatrix->setDimensionLabel(0,i,l_pInputMatrix->getDimensionLabel(0,i));
			l_pOutputPhaseMatrix->setDimensionLabel(0,i,l_pInputMatrix->getDimensionLabel(0,i));
		}

	}

	if(this->isInputTriggerActive(OVP_Algorithm_HilbertTransform_InputTriggerId_Process))
	{

		//Computing Hilbert transform for all channels
		for(uint32 channel=0; channel<l_ui32ChannelCount; channel++)
		{
			//Initialization of buffer vectors
			m_vecXcdSignalBuffer = VectorXcd::Zero(l_ui32SamplesPerChannel);
			m_vecXcdSignalFourier = VectorXcd::Zero(l_ui32SamplesPerChannel);

			//Initialization of vector h used to compute analytic signal
			m_vecXdHilbert.resize(l_ui32SamplesPerChannel);
			m_vecXdHilbert(0) = 1.0;

			if(l_ui32SamplesPerChannel%2 == 0)
			{
				m_vecXdHilbert(l_ui32SamplesPerChannel/2) = 1.0;

				m_vecXdHilbert.segment(1,(l_ui32SamplesPerChannel/2)-1).setOnes();
				m_vecXdHilbert.segment(1,(l_ui32SamplesPerChannel/2)-1) *= 2.0;

				m_vecXdHilbert.tail((l_ui32SamplesPerChannel/2)+1).setZero();
			}
			else
			{
				m_vecXdHilbert((l_ui32SamplesPerChannel+1)/2) = 1.0;

				m_vecXdHilbert.segment(1,(l_ui32SamplesPerChannel/2)).setOnes();
				m_vecXdHilbert.segment(1,(l_ui32SamplesPerChannel/2)) *= 2.0;

				m_vecXdHilbert.tail(((l_ui32SamplesPerChannel+1)/2)+1).setZero();
			}

			//Copy input signal chunk on buffer
			for(uint32 samples=0; samples<l_ui32SamplesPerChannel;samples++)
			{
				m_vecXcdSignalBuffer(samples).real(l_pInputMatrix->getBuffer()[samples + channel * (l_ui32SamplesPerChannel)]);
				m_vecXcdSignalBuffer(samples).imag(0.0);
			}

			//Fast Fourier Transform of input signal
			m_fft.fwd(m_vecXcdSignalFourier, m_vecXcdSignalBuffer);

			//Apply Hilbert transform by element-wise multiplying fft vector by h
			m_vecXcdSignalFourier = m_vecXcdSignalFourier.cwiseProduct(m_vecXdHilbert);


			//Inverse Fast Fourier transform
			m_fft.inv(m_vecXcdSignalBuffer, m_vecXcdSignalFourier); // m_vecXcdSignalBuffer is now the analytical signal of the initial input signal

			//Compute envelope and phase and pass it to the corresponding output
			for(uint32 samples=0; samples<l_ui32SamplesPerChannel;samples++)
			{

				l_pOutputHilbertMatrix->getBuffer()[samples + channel*l_ui32SamplesPerChannel] = m_vecXcdSignalBuffer(samples).imag();
				l_pOutputEnvelopeMatrix->getBuffer()[samples + channel*l_ui32SamplesPerChannel] = abs(m_vecXcdSignalBuffer(samples));
				l_pOutputPhaseMatrix->getBuffer()[samples + channel*l_ui32SamplesPerChannel] = arg(m_vecXcdSignalBuffer(samples));
			}

		}

	}
	return true;
}

#endif //TARGET_HAS_ThirdPartyEIGEN