implentation of CTCFakctor and WIC factor

include/ITANCTC.h

@@ -125,10 +125,16 @@ public:
 		class ITA_CTC_API Loudspeaker
 		{
 		public:
+			enum
+			{
+				LEFT_SIDE = 1,
+				RIGHT_SIDE = 2
+			};
 			Loudspeaker();
 			Loudspeaker( const Pose& oStaticPose );
 			Pose oPose;							//!< Pose of loudspeaker
 			const DAFFContentDFT* pDirectivity;	//!< Directivity data
+			int iSide;
 		};
 		std::vector< Loudspeaker > voLoudspeaker; //!< Loudspeaker of this setup
 
@@ -162,6 +168,7 @@ public:
 	/**
 	  * \return Static configuration parameters
 	  */
+
 	const Config& GetConfig() const;
 
 	//! Get beta parameter
@@ -309,7 +316,7 @@ private:
 	std::vector< ITAHDFTSpectra* > m_vpHelper2x2;	//!< Two-by-two helper matrix
 	ITAHDFTSpectrum* t;		//!< Helper
 	ITAHDFTSpectrum* det;	//!< Helper
-
+	int GetLoudspeakerSide(int);
 	//! Adds a HRIR into the target filter
 	/**
 	* @param oLoudspeakerPose	Pose of loudspeaker (or virtual source)

src/ITANCTC.cpp

@@ -79,10 +79,56 @@ int ITANCTC::GetN() const
 void ITANCTC::UpdateHeadPose( const Pose& oHead )
 {
     m_oHeadPose = oHead;
-
+	
 	return;
 }
 
+int ITANCTC::GetLoudspeakerSide(int iNumLoudspeaker)
+{
+	Pose oDest = GetLoudspeakerPose(iNumLoudspeaker);
+	VistaVector3D vConn = oDest.vPos - m_oHeadPose.vPos;
+	double dDistanceMeters = double(vConn.GetLength());
+	float fLS2HeadDelaySamples = (float)dDistanceMeters / m_oConfig.fSpeedOfSound * (float)m_oConfig.dSampleRate;
+
+	// 1/r attenuation due to distance law
+	// Note: Gain of HRIR dataset is normalized to 1 m according to convention
+	float fDistanceGain = 1 / (float)dDistanceMeters;
+
+	// @todo jst: mit VistaQuaternion lsen
+	/*
+	VistaVector3D oFrom = vConn;
+	oFrom.Normalize();
+	VistaQuaternion qOri( Vista::ViewVector, oFrom );
+	VistaQuaternion qHead2LS = qOri * m_oHeadPose.qOrient;
+	float fX, fY, fZ;
+	qHead2LS.GetAngles( fX, fY, fZ );
+
+	VistaVector3D v3View = qHead2LS.GetViewDir();
+
+	float fPhi = 180.0f / PI_F * fY;
+	float fTheta = 180.0f / PI_F * fX;
+	*/
+
+	double dPhiDeg, dThetaDeg;
+
+	VistaVector3D vDir = oDest.vPos - m_oHeadPose.vPos;
+	vDir.Normalize();
+	VistaVector3D vViewMinusZ = m_oHeadPose.vView * (-1.0f); // local z axis
+	const VistaVector3D vRight = vViewMinusZ.Cross(m_oHeadPose.vUp); // local x axis
+	const double dAzimuthAngleDeg = atan2(vDir.Dot(vRight), vDir.Dot(m_oHeadPose.vView)) * 180.0f / ITAConstants::PI_D;
+	dPhiDeg = ((dAzimuthAngleDeg < 0.0f) ? (dAzimuthAngleDeg + 360.0f) : dAzimuthAngleDeg);
+	dThetaDeg = asin(vDir.Dot(m_oHeadPose.vUp)) * 180.0f / ITAConstants::PI_D;
+	if (dPhiDeg < 180)
+	{
+		return Config::Loudspeaker::LEFT_SIDE;
+	}
+	else
+	{
+		return Config::Loudspeaker::RIGHT_SIDE;
+	}
+	
+}
+
 bool ITANCTC::AddHRIR( const Pose& oDest, ITASampleFrame& sfDestHRIR, bool& bOutOfRange, double dReflectionFactor/* = 1.0f */ ) const
 {
 	if( sfDestHRIR.channels() != 2 )
@@ -121,7 +167,7 @@ bool ITANCTC::AddHRIR( const Pose& oDest, ITASampleFrame& sfDestHRIR, bool& bOut
 	const double dAzimuthAngleDeg = atan2( vDir.Dot( vRight ), vDir.Dot( m_oHeadPose.vView ) ) * 180.0f / ITAConstants::PI_D;
 	dPhiDeg = ( ( dAzimuthAngleDeg < 0.0f ) ? ( dAzimuthAngleDeg + 360.0f ) : dAzimuthAngleDeg );
 	dThetaDeg = asin( vDir.Dot( m_oHeadPose.vUp ) ) * 180.0f / ITAConstants::PI_D;
-
+	
 	int iHRIRPreOffset = m_pHRIR->getMinEffectiveFilterOffset();
 	int iHRIRFilerTaps = m_pHRIR->getMaxEffectiveFilterLength();
 	int iOffset = (std::max)( 0, int( fLS2HeadDelaySamples - iHRIRPreOffset ) ); // remove starting zeros from HRIR
@@ -198,27 +244,60 @@ bool ITANCTC::CalculateFilter( std::vector< ITAHDFTSpectra* >& vpCTCFilter )
 	for( int n=0; n<GetN(); n++ )
 	{
 		ITAHDFTSpectra* pHRTF( m_vpHRTFs[n] ); // two-channel
+
 		float fWeight = float( m_vdWeights[n] ); // diag element
+		
+		float fLeftChannel = 1;
+		float fRightChannel = 1;
+
+		float fLeftChannelWICK = 1;
+		float fRightChannelWICK = 1;
 
+		if (GetLoudspeakerSide(n) == Config::Loudspeaker::LEFT_SIDE)
+		{
+			fRightChannel *= m_fCrossTalkCancellationFactor;
+			fLeftChannelWICK *= m_fWaveIncidenceAngleCompensationFactor;
+		}
+		else if (GetLoudspeakerSide(n) == Config::Loudspeaker::RIGHT_SIDE)
+		{
+			fLeftChannel *= m_fCrossTalkCancellationFactor;
+			fRightChannelWICK *= m_fWaveIncidenceAngleCompensationFactor;
+		}
 		// Element wise (a and d): HWH* -> 2x2
 		t->copy( (*pHRTF)[0] );
-		t->mul( fWeight );
+		t->mul( fWeight*fLeftChannel );
+		t->log();
+		t->mul(fLeftChannelWICK);
+		t->exp();
+		
 		t->mul_conj( (*pHRTF)[0] );
 		n == 0 ? a->copy( t ) : a->add( t );
 
 		t->copy( (*pHRTF)[1] );
-		t->mul( fWeight );
+		t->mul( fWeight*fRightChannel );
+		t->log();
+		t->mul(fRightChannelWICK);
+		t->exp();
+
 		t->mul_conj( (*pHRTF)[1] );
 		n == 0 ? d->copy( t ) : d->add( t );
-
+		
 		// Cross elements (b and c): HWH*
 		t->copy( (*pHRTF)[1] );		
-		t->mul( fWeight );
+		t->mul( fWeight*fRightChannel );
+		t->log();
+		t->mul(fRightChannelWICK);
+		t->exp();
+
 		t->mul_conj( (*pHRTF)[0] );
 		n == 0 ? b->copy( t ) : b->add( t );
 
 		t->copy( (*pHRTF)[0] );
-		t->mul( fWeight );
+		t->mul( fWeight*fLeftChannel );
+		t->log();
+		t->mul(fLeftChannelWICK);
+		t->exp();
+
 		t->mul_conj( (*pHRTF)[1] );
 		n == 0 ? c->copy( t ) : c->add( t );
 	}
@@ -260,26 +339,44 @@ bool ITANCTC::CalculateFilter( std::vector< ITAHDFTSpectra* >& vpCTCFilter )
 
 	for( int n=0; n<GetN(); n++ )
 	{
+		float fLeftChannel = 1;
+		float fRightChannel = 1;
+
+		float fLeftChannelWICK = 1;
+		float fRightChannelWICK = 1;
+
+		if (GetLoudspeakerSide(n) == Config::Loudspeaker::LEFT_SIDE)
+		{
+			fRightChannel *= m_fCrossTalkCancellationFactor;
+			fLeftChannelWICK *= m_fWaveIncidenceAngleCompensationFactor;
+		}
+		else if (GetLoudspeakerSide(n) == Config::Loudspeaker::RIGHT_SIDE)
+		{
+			fLeftChannel *= m_fCrossTalkCancellationFactor;
+			fRightChannelWICK *= m_fWaveIncidenceAngleCompensationFactor;
+		}
+
 		ITAHDFTSpectra* pHRTF( m_vpHRTFs[n] ); // two-channel
 		ITAHDFTSpectra* pCTCFilter( vpCTCFilter[n] ); // two-channel
 		float fWeight = float( m_vdWeights[n] ); // diag element
 
 		t->copy( a );
 		t->mul_conj( (*pHRTF)[0] );
-		t->mul( fWeight );
+		t->mul( fWeight*fLeftChannel );
 		(*pCTCFilter)[0]->copy( t );
+
 		t->copy( b );
 		t->mul_conj( (*pHRTF)[1] );
-		t->mul( fWeight );
+		t->mul( fWeight*fRightChannel );
 		(*pCTCFilter)[0]->add( t );
 
 		t->copy( c );
 		t->mul_conj( (*pHRTF)[0] );
-		t->mul( fWeight );
+		t->mul( fWeight*fLeftChannel );
 		(*pCTCFilter)[1]->copy( t );
 		t->copy( d );
 		t->mul_conj( (*pHRTF)[1] );
-		t->mul( fWeight );
+		t->mul( fWeight*fRightChannel );
 		(*pCTCFilter)[1]->add( t );
 
 		//pCTCFilter->Export( "CTCFilter_noshift_" + IntToString( n+1 ) );
@@ -433,12 +530,14 @@ ITANCTC::Pose ITANCTC::GetHeadPose() const
 // --- Loudspeaker ---
 
 ITANCTC::Config::Loudspeaker::Loudspeaker()
-	: pDirectivity( NULL )
+: pDirectivity(NULL)
+, iSide(0)
 {
 }
 
 ITANCTC::Config::Loudspeaker::Loudspeaker( const Pose& oStaticPose )
 	: pDirectivity( NULL )
+	, iSide(0)
 {
 	oPose = oStaticPose;
 }