ITANCTC.cpp 17.4 KB
Newer Older
Jonas Stienen's avatar
Jonas Stienen committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
// $Id: ITANCTC.cpp 2395 2012-04-20 06:58:52Z stienen $

#include <ITANCTC.h>
#include <ITACTCUtils.h>

#include <assert.h>
#include <complex>
#include <cmath>
#include <iostream>

#include <DAFF.h>

#include <ITAConstants.h>
#include <ITAException.h>
#include <ITAHDFTSpectrum.h>
#include <ITANumericUtils.h>
#include <ITAStringUtils.h>

ITANCTC::ITANCTC( const Config& oNCTCConfig )
	: m_oConfig( oNCTCConfig )
	, m_pHRIR( NULL )
{
	m_fBeta = float( 1e-4 );
24

25
26
27
	if( oNCTCConfig.N <= 0 )
		ITA_EXCEPT1( INVALID_PARAMETER, "Wrong configuration, N-CTC requires a valid number of loudspeakers" );

Jonas Stienen's avatar
Jonas Stienen committed
28
	m_iOptimization = m_oConfig.iOptimization;
29
30
	m_fCrossTalkCancellationFactor = m_oConfig.fCrossTalkCancellationFactor;
	m_fWaveIncidenceAngleCompensationFactor = m_oConfig.fWaveIncidenceAngleCompensationFactor;
Jonas Stienen's avatar
Jonas Stienen committed
31
32
33

	m_oHeadPose.vPos.SetToZeroVector();
	m_oHeadPose.vView.SetValues( 0, 0, -1.0f );
34
	m_oHeadPose.vUp.SetValues( 0, 1.0f, 0 );
Jonas Stienen's avatar
Jonas Stienen committed
35
36
37
	//m_oHeadPose.qOrient.SetToNeutralQuaternion();

	m_sfCTC_temp.init( 2, m_oConfig.iCTCFilterLength, true );
38
39
40

	int iDFTSize = m_oConfig.iCTCFilterLength + 1;
	for( int n = 0; n < GetN(); n++ )
Jonas Stienen's avatar
Jonas Stienen committed
41
42
43
	{
		m_vdWeights.push_back( 1.0f );
		m_vpHRTFs.push_back( new ITAHDFTSpectra( m_oConfig.dSampleRate, 2, iDFTSize, true ) );
44
		m_vfDelayTime.push_back( float( m_oConfig.iCTCFilterLength / m_oConfig.dSampleRate / 2.0f ) );
Jonas Stienen's avatar
Jonas Stienen committed
45
46
	}

47
	for( int i = 0; i < 2; i++ )
Jonas Stienen's avatar
Jonas Stienen committed
48
49
50
51
52
		m_vpHelper2x2.push_back( new ITAHDFTSpectra( m_oConfig.dSampleRate, 2, iDFTSize, true ) );

	t = new ITAHDFTSpectrum( m_oConfig.dSampleRate, iDFTSize, true );
	det = new ITAHDFTSpectrum( m_oConfig.dSampleRate, iDFTSize, true );

53
54
55
	int l = m_oConfig.iCTCFilterLength;
	m_fft.plan( ITAFFT::FFT_R2C, l, m_sfCTC_temp[ 0 ].data(), ( *m_vpHRTFs[ 0 ] )[ 0 ]->data() );
	m_ifft.plan( ITAFFT::IFFT_C2R, l, ( *m_vpHRTFs[ 0 ] )[ 0 ]->data(), m_sfCTC_temp[ 0 ].data() );
Jonas Stienen's avatar
Jonas Stienen committed
56
57
58
59
60

}

ITANCTC::~ITANCTC()
{
61
62
	for( int n = 0; n < GetN(); n++ )
		delete m_vpHRTFs[ n ];
Jonas Stienen's avatar
Jonas Stienen committed
63

64
65
	for( int i = 0; i < 2; i++ )
		delete m_vpHelper2x2[ i ];
Jonas Stienen's avatar
Jonas Stienen committed
66
67

	delete t, det;
68

69
	return;
Jonas Stienen's avatar
Jonas Stienen committed
70
71
72
73
74
75
76
77
78
79
80
81
82
83
}

const ITANCTC::Config& ITANCTC::GetConfig() const
{
	return m_oConfig;
}

int ITANCTC::GetN() const
{
	return m_oConfig.N;
}

void ITANCTC::UpdateHeadPose( const Pose& oHead )
{
84
85
	m_oHeadPose = oHead;

Jonas Stienen's avatar
Jonas Stienen committed
86
87
88
	return;
}

89
int ITANCTC::GetLoudspeakerSide( int iNumLoudspeaker )
90
{
91
	Pose oDest = GetLoudspeakerPose( iNumLoudspeaker );
92
	VistaVector3D vConn = oDest.vPos - m_oHeadPose.vPos;
93
94
	double dDistanceMeters = double( vConn.GetLength() );
	float fLS2HeadDelaySamples = ( float ) dDistanceMeters / m_oConfig.fSpeedOfSound * ( float ) m_oConfig.dSampleRate;
95
96
97

	// 1/r attenuation due to distance law
	// Note: Gain of HRIR dataset is normalized to 1 m according to convention
98
	float fDistanceGain = 1 / ( float ) dDistanceMeters;
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118

	// @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();
119
120
121
122
123
124
	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 )
125
126
127
128
129
130
131
	{
		return Config::Loudspeaker::LEFT_SIDE;
	}
	else
	{
		return Config::Loudspeaker::RIGHT_SIDE;
	}
132

133
134
}

Jonas Stienen's avatar
Jonas Stienen committed
135
136
137
138
139
140
141
142
bool ITANCTC::AddHRIR( const Pose& oDest, ITASampleFrame& sfDestHRIR, bool& bOutOfRange, double dReflectionFactor/* = 1.0f */ ) const
{
	if( sfDestHRIR.channels() != 2 )
		ITA_EXCEPT1( INVALID_PARAMETER, "Two channel HRIR expected" );

	// Calculate sample delay of HRIR insertion position
	VistaVector3D vConn = oDest.vPos - m_oHeadPose.vPos;
	double dDistanceMeters = double( vConn.GetLength() );
143
144
	float fLS2HeadDelaySamples = ( float ) dDistanceMeters / m_oConfig.fSpeedOfSound * ( float ) m_oConfig.dSampleRate;

Jonas Stienen's avatar
Jonas Stienen committed
145
146
	// 1/r attenuation due to distance law
	// Note: Gain of HRIR dataset is normalized to 1 m according to convention
147
	float fDistanceGain = 1 / ( float ) dDistanceMeters;
Jonas Stienen's avatar
Jonas Stienen committed
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167

	// @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();
168
	VistaVector3D vViewMinusZ = m_oHeadPose.vView * ( -1.0f ); // local z axis
Jonas Stienen's avatar
Jonas Stienen committed
169
170
171
172
	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;
173

Jonas Stienen's avatar
Jonas Stienen committed
174
175
	int iHRIRPreOffset = m_pHRIR->getMinEffectiveFilterOffset();
	int iHRIRFilerTaps = m_pHRIR->getMaxEffectiveFilterLength();
176
177
	int iOffset = ( std::max )( 0, int( fLS2HeadDelaySamples - iHRIRPreOffset ) ); // remove starting zeros from HRIR

Jonas Stienen's avatar
Jonas Stienen committed
178
179
180
181
182
183
184
185
	iOffset = int( fLS2HeadDelaySamples ) + iHRIRPreOffset;

	// Check against buffer overrun, prevent if necessary (PROBLEM: LS too far away for target filter)
	if( iOffset + iHRIRFilerTaps > sfDestHRIR.length() )
		return false;

	int iRecordIndex;
	m_pHRIR->getNearestNeighbour( DAFF_OBJECT_VIEW, float( dPhiDeg ), float( dThetaDeg ), iRecordIndex, bOutOfRange );
186
187
188
	m_pHRIR->addFilterCoeffs( iRecordIndex, 0, sfDestHRIR[ 0 ].data() + iOffset, fDistanceGain * float( dReflectionFactor ) );
	m_pHRIR->addFilterCoeffs( iRecordIndex, 1, sfDestHRIR[ 1 ].data() + iOffset, fDistanceGain * float( dReflectionFactor ) );

Jonas Stienen's avatar
Jonas Stienen committed
189
190
191
192
193
194
195
196
	return true;
}

const ITANCTC::Pose& ITANCTC::GetLoudspeakerPose( int iLoudspeakerID ) const
{
	if( iLoudspeakerID > GetN() )
		ITA_EXCEPT1( INVALID_PARAMETER, "Loudspeaker ID (starting from 1) out of range." );

197
	return m_oConfig.voLoudspeaker[ iLoudspeakerID - 1 ].oPose;
Jonas Stienen's avatar
Jonas Stienen committed
198
199
200
201
202
203
}

bool ITANCTC::CalculateFilter( std::vector< ITAHDFTSpectra* >& vpCTCFilter )
{
	if( !m_pHRIR )
		return false;
204

Jonas Stienen's avatar
Jonas Stienen committed
205
	bool bOutOfRange = false;
206
	for( int n = 0; n < GetN(); n++ )
Jonas Stienen's avatar
Jonas Stienen committed
207
208
209
210
211
	{
		m_sfCTC_temp.zero();

		switch( m_iOptimization )
		{
212
			// Here comes more if ready ... see feature/room_compensation branch
213
		case Config::OPTIMIZATION_NONE:
214
215
216
217
		default:
			if( AddHRIR( GetLoudspeakerPose( n + 1 ), m_sfCTC_temp, bOutOfRange ) == false )
				return false;
			break;
Jonas Stienen's avatar
Jonas Stienen committed
218
219
220
221
222
		}

		if( bOutOfRange )
			return false;

223
		ITAHDFTSpectra* pHRTF( m_vpHRTFs[ n ] );
Jonas Stienen's avatar
Jonas Stienen committed
224
225

		// Convert HRIRs to HRTFs
226
227
		m_fft.execute( m_sfCTC_temp[ 0 ].data(), ( *pHRTF )[ 0 ]->data() );
		m_fft.execute( m_sfCTC_temp[ 1 ].data(), ( *pHRTF )[ 1 ]->data() );
Jonas Stienen's avatar
Jonas Stienen committed
228

229
230
231
#ifdef NCTC_EXPORT_FILTER_TO_HARDDRIVE
		pHRTF->Export( "HRIR_LS" + IntToString( n+1 ) + "_RAW" );
#endif // NCTC_EXPORT_FILTER_TO_HARDDRIVE
Jonas Stienen's avatar
Jonas Stienen committed
232
	}
233

Jonas Stienen's avatar
Jonas Stienen committed
234
235
236
237
	/* Matrix to be inverted (as two-by-two row vector)
	 * a  c
	 * b  d
	 */
238
239
240
241
	ITAHDFTSpectrum* a = ( *m_vpHelper2x2[ 0 ] )[ 0 ];
	ITAHDFTSpectrum* b = ( *m_vpHelper2x2[ 0 ] )[ 1 ];
	ITAHDFTSpectrum* c = ( *m_vpHelper2x2[ 1 ] )[ 0 ];
	ITAHDFTSpectrum* d = ( *m_vpHelper2x2[ 1 ] )[ 1 ];
Jonas Stienen's avatar
Jonas Stienen committed
242
243

	// Least-squares minimization: C = WH*(HWH*-\beta)^-1
244
245
	// using H* as the hermitian (complex conjugated) transpose of H

246
247
248
249
250
251
252
	/* Sweet spot optimization: modify input HRTFs
	 *   1. Crosstalk side will be lowered by CTC factor
	 *   2. HRTF will be flattened by WICK factor
	 */ 


	// Prepare HRTF spectra
Jonas Stienen's avatar
Jonas Stienen committed
253

254
255
	int iDFTSize = m_oConfig.iCTCFilterLength + 1;
	for( int n = 0; n < GetN(); n++ )
Jonas Stienen's avatar
Jonas Stienen committed
256
	{
257
258
		ITAHDFTSpectra* pHRTF( m_vpHRTFs[ n ] ); // two-channel

259

260
261
262
263
264
265
		//  --- WICK factor ---
		
		// First, store original energy of HRTF (left channel)
		float fEnergy = ( *pHRTF )[ 0 ]->getEnergy();
		
		// Apply WICK factor in cepstrum domain (left channel)
266
267
268
269
		( *pHRTF )[ 0 ]->log();
		( *pHRTF )[ 0 ]->mul( m_fWaveIncidenceAngleCompensationFactor );
		( *pHRTF )[ 0 ]->exp();

270
		// Compensate initial HRTF energy when WICK is used (left channel)
271
272
273
274
275
		assert( fEnergy > 0 );
		float fEnergyCompensation = powf( fEnergy, ( 1 - m_fWaveIncidenceAngleCompensationFactor ) );
		( *pHRTF )[ 0 ]->mul( fEnergyCompensation );


276
		// First, store original energy of HRTF (right channel)
277
		fEnergy = ( *pHRTF )[ 1 ]->getEnergy();
278
279
280
281
282
283
284

		// Apply WICK factor in cepstrum domain (right channel)
		( *pHRTF )[ 1 ]->log();
		( *pHRTF )[ 1 ]->mul( m_fWaveIncidenceAngleCompensationFactor );
		( *pHRTF )[ 1 ]->exp();

		// Compensate initial HRTF energy when WICK is used (right channel)
285
286
287
		assert( fEnergy > 0 );
		fEnergyCompensation = powf( fEnergy, ( 1 - m_fWaveIncidenceAngleCompensationFactor ) );
		( *pHRTF )[ 1 ]->mul( fEnergyCompensation );
288
		
289

290
		//  --- CTC compensation factor + weighting ---
291
		
292
293
		// Weighting
		float fWeight = float( m_vdWeights[ n ] ); // diag element
294

295
296
297
298
299
300
301
302
303
304
		// CTC compensation factors
		int iLSSide = GetLoudspeakerSide( n + 1 );

		float fRightChannelCTC = 1.0f;
		if( iLSSide == Config::Loudspeaker::LEFT_SIDE )
			fRightChannelCTC *= m_fCrossTalkCancellationFactor; // only apply to left channel if LS is at left side

		float fLeftChannelCTC = 1.0f;
		if( iLSSide == Config::Loudspeaker::RIGHT_SIDE )
			fLeftChannelCTC *= m_fCrossTalkCancellationFactor; // only apply to right channel if LS is at right side
Jonas Stienen's avatar
Jonas Stienen committed
305
306

		// Element wise (a and d): HWH* -> 2x2
307
308
309
310
		t->copy( ( *pHRTF )[ 0 ] );
		t->mul( fWeight * fLeftChannelCTC );

		t->mul_conj( ( *pHRTF )[ 0 ] );
Jonas Stienen's avatar
Jonas Stienen committed
311
312
		n == 0 ? a->copy( t ) : a->add( t );

313
314
315
316
		t->copy( ( *pHRTF )[ 1 ] );
		t->mul( fWeight * fRightChannelCTC );

		t->mul_conj( ( *pHRTF )[ 1 ] );
Jonas Stienen's avatar
Jonas Stienen committed
317
		n == 0 ? d->copy( t ) : d->add( t );
318

Jonas Stienen's avatar
Jonas Stienen committed
319
		// Cross elements (b and c): HWH*
320
321
		t->copy( ( *pHRTF )[ 1 ] );
		t->mul( fWeight * fRightChannelCTC );
322

323
		t->mul_conj( ( *pHRTF )[ 0 ] );
Jonas Stienen's avatar
Jonas Stienen committed
324
325
		n == 0 ? b->copy( t ) : b->add( t );

326
327
328
329
		t->copy( ( *pHRTF )[ 0 ] );
		t->mul( fWeight * fLeftChannelCTC );

		t->mul_conj( ( *pHRTF )[ 1 ] );
Jonas Stienen's avatar
Jonas Stienen committed
330
331
332
		n == 0 ? c->copy( t ) : c->add( t );
	}

333

Jonas Stienen's avatar
Jonas Stienen committed
334
	// Regularize
335

Jonas Stienen's avatar
Jonas Stienen committed
336
337
338
339
	a->add( m_fBeta );
	d->add( m_fBeta );

	ITAHDFTSpectra abcd( m_oConfig.dSampleRate, 4, m_oConfig.iCTCFilterLength );
340
341
342
343
	abcd[ 0 ]->copyFrom( *a );
	abcd[ 1 ]->copyFrom( *b );
	abcd[ 2 ]->copyFrom( *c );
	abcd[ 3 ]->copyFrom( *d );
Jonas Stienen's avatar
Jonas Stienen committed
344

345
346

	// Invert via determinant (simple 2x2 matrix inversion)
Jonas Stienen's avatar
Jonas Stienen committed
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
	/*
	 * d/det  -b/det
	 * -c/det  a/det
	 */
	det->copy( a );
	det->mul( d );
	t->copy( b );
	t->mul( c );
	det->sub( t );

	t->copy( d );
	d->copy( a );
	a->copy( t );
	b->mul( -1.0f );
	c->mul( -1.0f );
	a->div( det );
	b->div( det );
	c->div( det );
	d->div( det );


	// Calculate CTC filter WH*Inv (Nx2)
369

Jonas Stienen's avatar
Jonas Stienen committed
370
371
	ITASampleFrame sfTargetData_shift( m_sfCTC_temp.channels(), m_sfCTC_temp.length(), true );

372
	for( int n = 0; n < GetN(); n++ )
Jonas Stienen's avatar
Jonas Stienen committed
373
	{
374
375
		// Again, apply CTC compensation factors
		int iLSSide = GetLoudspeakerSide( n + 1 );
376

377
378
379
		float fRightChannelCTC = 1.0f;
		if( iLSSide == Config::Loudspeaker::LEFT_SIDE )
			fRightChannelCTC *= m_fCrossTalkCancellationFactor;
380

381
382
383
		float fLeftChannelCTC = 1.0f;
		if( iLSSide == Config::Loudspeaker::RIGHT_SIDE )
			fLeftChannelCTC *= m_fCrossTalkCancellationFactor;
384

385
386
387
		ITAHDFTSpectra* pHRTF( m_vpHRTFs[ n ] ); // two-channel, already WICKed
		ITAHDFTSpectra* pCTCFilter( vpCTCFilter[ n ] ); // two-channel
		float fWeight = float( m_vdWeights[ n ] ); // diag element
Jonas Stienen's avatar
Jonas Stienen committed
388
389

		t->copy( a );
390
391
392
		t->mul_conj( ( *pHRTF )[ 0 ] );
		t->mul( fWeight * fLeftChannelCTC );
		( *pCTCFilter )[ 0 ]->copy( t );
393

Jonas Stienen's avatar
Jonas Stienen committed
394
		t->copy( b );
395
396
397
		t->mul_conj( ( *pHRTF )[ 1 ] );
		t->mul( fWeight * fRightChannelCTC );
		( *pCTCFilter )[ 0 ]->add( t );
Jonas Stienen's avatar
Jonas Stienen committed
398
399

		t->copy( c );
400
401
402
		t->mul_conj( ( *pHRTF )[ 0 ] );
		t->mul( fWeight * fLeftChannelCTC );
		( *pCTCFilter )[ 1 ]->copy( t );
Jonas Stienen's avatar
Jonas Stienen committed
403
		t->copy( d );
404
405
406
		t->mul_conj( ( *pHRTF )[ 1 ] );
		t->mul( fWeight * fRightChannelCTC );
		( *pCTCFilter )[ 1 ]->add( t );
Jonas Stienen's avatar
Jonas Stienen committed
407

408
409
410
411
#ifdef NCTC_EXPORT_FILTER_TO_HARDDRIVE
		pCTCFilter->Export( "CTCFilter_noshift_" + IntToString( n+1 ) );
#endif // NCTC_EXPORT_FILTER_TO_HARDDRIVE 

Jonas Stienen's avatar
Jonas Stienen committed
412
413

		// Time-shift
414
415
416

		m_ifft.execute( ( *pCTCFilter )[ 0 ]->data(), m_sfCTC_temp[ 0 ].data() );
		m_ifft.execute( ( *pCTCFilter )[ 1 ]->data(), m_sfCTC_temp[ 1 ].data() );
Jonas Stienen's avatar
Jonas Stienen committed
417
418
419
420
421

		// Normalize after IFFT
		m_sfCTC_temp.div_scalar( float( m_sfCTC_temp.length() ) );

		// Shift the CTC filter according to desired delay
422
423
		int iShiftSamples = int( m_oConfig.dSampleRate * m_vfDelayTime[ n ] );
		if( m_vfDelayTime[ n ] < 0.0f )
Jonas Stienen's avatar
Jonas Stienen committed
424
425
			iShiftSamples = m_sfCTC_temp.length() / 2; // if invalid, shift by half length
		sfTargetData_shift.cyclic_write( m_sfCTC_temp, m_sfCTC_temp.length(), 0, iShiftSamples );
426

427
428
		m_fft.execute( sfTargetData_shift[ 0 ].data(), ( *pCTCFilter )[ 0 ]->data() );
		m_fft.execute( sfTargetData_shift[ 1 ].data(), ( *pCTCFilter )[ 1 ]->data() );
Jonas Stienen's avatar
Jonas Stienen committed
429

430
431
432
#ifdef NCTC_EXPORT_FILTER_TO_HARDDRIVE
		pCTCFilter->Export( "CTCFilter_shift_" + IntToString( n+1 ) );
#endif // NCTC_EXPORT_FILTER_TO_HARDDRIVE
Jonas Stienen's avatar
Jonas Stienen committed
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
	}

	return true;
}

void ITANCTC::SetHRIR( const DAFFContentIR* pHRIR )
{
	// TODO: Set m_oConfig.iCTCFilterLength according to HRIR length and working area
	// TODO: Minimize loudspeaker delay

	if( m_pHRIR != pHRIR )
	{
		m_pHRIR = pHRIR;

		// Return if HRIR is reset (NULL)
448
		if( !m_pHRIR )
Jonas Stienen's avatar
Jonas Stienen committed
449
			return;
450
451
452
453

		if( m_pHRIR->getProperties()->getNumberOfChannels() != 2 )
			ITA_EXCEPT1( INVALID_PARAMETER, "HRIR dataset must have exactly two channels" );

Jonas Stienen's avatar
Jonas Stienen committed
454
455
		// Get HRIR filter length
		if( m_pHRIR->getFilterLength() > m_oConfig.iCTCFilterLength )
456
457
458
			ITA_EXCEPT1( INVALID_PARAMETER,
			std::string( "The filter length of the HRIR database \"" ) + m_pHRIR->getParent()->getFilename() +
			std::string( "\" is larger than the CTC filters. This leads to buffer overruns." ) );
Jonas Stienen's avatar
Jonas Stienen committed
459
460
461
462

		// Get HRIR delay samples
		const DAFFMetadata* pHRIRMetadata = m_pHRIR->getParent()->getMetadata();

463
464
465
466
467
468
469
470
471
472
473
		if( !pHRIRMetadata->hasKey( "DELAY_SAMPLES" ) )
			ITA_EXCEPT1( INVALID_PARAMETER,
			std::string( "HRIR database \"" ) + m_pHRIR->getParent()->getFilename() +
			std::string( "\" is missing a \"DELAY_SAMPLES\" metadata tag" ) );
		if( ( pHRIRMetadata->getKeyType( "DELAY_SAMPLES" ) != DAFFMetadata::DAFF_INT ) &&
			( pHRIRMetadata->getKeyType( "DELAY_SAMPLES" ) != DAFFMetadata::DAFF_FLOAT ) )
			ITA_EXCEPT1( INVALID_PARAMETER,
			std::string( "The metadata tag \"DELAY_SAMPLES\" in HRIR database \"" ) + m_pHRIR->getParent()->getFilename() +
			std::string( "\" is not numeric" ) );

		float fHRIRDelay = ( float ) pHRIRMetadata->getKeyFloat( "DELAY_SAMPLES" );
Jonas Stienen's avatar
Jonas Stienen committed
474
		if( fHRIRDelay < 0 )
475
476
477
			ITA_EXCEPT1( INVALID_PARAMETER,
			std::string( "The metadata tag \"DELAY_SAMPLES\" in HRIR database \"" ) + m_pHRIR->getParent()->getFilename() +
			std::string( "\" must not be negative" ) );
Jonas Stienen's avatar
Jonas Stienen committed
478
479
480
481
482
483
484
485
486


	}
}

void ITANCTC::SetBeta( float fBeta )
{
	m_fBeta = fBeta;
}
487

488
489
490
491
492
float ITANCTC::GetBeta()
{
	return m_fBeta;
}

493
void ITANCTC::SetCrossTalkCancellationFactor( float fFactor )
494
{
495
	m_fCrossTalkCancellationFactor = fFactor;
496
497
}

498
float ITANCTC::GetCrossTalkCancellationFactor()
499
{
500
501
502
503
504
505
506
507
508
509
510
	return m_fCrossTalkCancellationFactor;
}

void ITANCTC::SetWaveIncidenceAngleCompensationFactor( float fFactor )
{
	m_fWaveIncidenceAngleCompensationFactor = fFactor;
}

float ITANCTC::GetWaveIncidenceAngleCompensation()
{
	return m_fWaveIncidenceAngleCompensationFactor;
511
512
}

513
void ITANCTC::SetDelayTime( float fDelayTime )
Jonas Stienen's avatar
Jonas Stienen committed
514
{
515
516
517
518
519
520
521
522
523
524
525
	for( int n = 0; n < GetN(); n++ )
		m_vfDelayTime[ n ] = fDelayTime;
}

void ITANCTC::SetDelayTime( std::vector< float > vfDelayTime )
{
	if( int( vfDelayTime.size() ) != GetN() )
		ITA_EXCEPT1( INVALID_PARAMETER, "Provide as many delay values as channels for NCTC" );

	for( int n = 0; n < GetN(); n++ )
		m_vfDelayTime[ n ] = vfDelayTime[ n ];
Jonas Stienen's avatar
Jonas Stienen committed
526
527
}

528
std::vector< float > ITANCTC::GetDelayTime()
Jonas Stienen's avatar
Jonas Stienen committed
529
{
530
531
532
533
	std::vector< float > vfDelayTime( GetN() );
	for( int n = 0; n < GetN(); n++ )
		vfDelayTime[ n ] = m_vfDelayTime[ n ];
	return vfDelayTime;
Jonas Stienen's avatar
Jonas Stienen committed
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
}

void ITANCTC::SetOptimization( int iOptimization )
{
	m_iOptimization = iOptimization;
}

int ITANCTC::GetOptimization() const
{
	return  m_iOptimization;
}

bool ITANCTC::GetHRTF( std::vector< ITAHDFTSpectra* >& vpHRTF ) const
{
	if( m_vpHRTFs.empty() )
		return false;

	vpHRTF = m_vpHRTFs;
	return true;
}

ITANCTC::Pose ITANCTC::GetHeadPose() const
{
	return m_oHeadPose;
}

560

Jonas Stienen's avatar
Jonas Stienen committed
561
562
563
// --- Loudspeaker ---

ITANCTC::Config::Loudspeaker::Loudspeaker()
564
565
	: pDirectivity( NULL )
	, iSide( SIDE_NOT_SPECIFIED )
Jonas Stienen's avatar
Jonas Stienen committed
566
567
568
569
570
{
}

ITANCTC::Config::Loudspeaker::Loudspeaker( const Pose& oStaticPose )
	: pDirectivity( NULL )
571
	, iSide( SIDE_NOT_SPECIFIED )
Jonas Stienen's avatar
Jonas Stienen committed
572
573
574
575
576
{
	oPose = oStaticPose;
}


577
// --- Pose ---
Jonas Stienen's avatar
Jonas Stienen committed
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594

ITANCTC::Pose& ITANCTC::Pose::operator=( const ITANCTC::Pose& oPoseRHS )
{
	vPos = oPoseRHS.vPos;
	//qOrient = oPoseRHS.qOrient;
	vView = oPoseRHS.vView;
	vUp = oPoseRHS.vUp;

	return *this;
}

void ITANCTC::Pose::SetOrientationYPRdeg( double fYaw, double fPitch, double fRoll )
{
	double yaw = grad2rad( fYaw );
	double pitch = grad2rad( fPitch );
	double roll = grad2rad( fRoll );

595
596
597
	double sy = sin( yaw ), cy = cos( yaw );
	double sp = sin( pitch ), cp = cos( pitch );
	double sr = sin( roll ), cr = cos( roll );
Jonas Stienen's avatar
Jonas Stienen committed
598
599
600
601
602
603

	vView.SetValues( -sy*cp, sp, -cy*cp );
	vUp.SetValues( cy*sr + sy*sp*cr, cp*cr, -sy*sr + cy*sp*cr );

	return;
}
604
605
606
607
608
609
610
611
612
613
614

ITANCTC::Config::Config()
{
	// Set some default values
	N = 0;
	iCTCFilterLength = 4096;
	fSpeedOfSound = 344.0f;
	iOptimization = OPTIMIZATION_NONE;
	fCrossTalkCancellationFactor = 1.0f;
	fWaveIncidenceAngleCompensationFactor = 1.0f;
}