ITANumericUtils.h 15.6 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
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
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
/*
* ----------------------------------------------------------------
*
*		ITA core libs
*		(c) Copyright Institute of Technical Acoustics (ITA)
*		RWTH Aachen University, Germany, 2015-2016
*
* ----------------------------------------------------------------
*				    ____  __________  _______
*				   //  / //__   ___/ //  _   |
*				  //  /    //  /    //  /_|  |
*				 //  /    //  /    //  ___   |
*				//__/    //__/    //__/   |__|
*
* ----------------------------------------------------------------
*
*/
// $Id: ITAStringUtils.h 3923 2015-01-07 16:09:21Z pelzer $

#ifndef INCLUDE_WATCHER_ITA_NUMERIC_UTILS
#define INCLUDE_WATCHER_ITA_NUMERIC_UTILS

#include <ITABaseDefinitions.h>

#include <algorithm>
#include <cmath>
#include <limits>
#include <vector>

// Ein Makro für die Berechnung von Prozentsätzen (Wieviel Prozent p sind P von G?)
// Wichtiges Feature: Behandlung des Falles G = 0. Rückgabewert in diesem Fall 0.
#define PERCENT_OF(P, G) (G == 0 ? 0 : 100*(double) P/(double) G)

//! Zweierexponent einer natürlichen Zahl ermitteln
/**
 * Ermittelt die natürliche Zahl k, welche 2^k = x erfüllt.
 * Existiert keine solche, d.h. ist x keine Zweierpotenz, so
 * gibt die Funktion -1 zurück.
 */
ITA_BASE_API int getExp2(unsigned int x);

//! Test auf Zweierpotenz
ITA_BASE_API bool isPow2(unsigned int x);

//! Nächstgrößeres Zweierpotenz zurückgeben
ITA_BASE_API unsigned int nextPow2(unsigned int x);

// Rundung auf Compilern definieren die dies nicht standardmäßig mitliefern

#ifndef HAVE_ROUND
ITA_BASE_API inline double round(double x) { return (x<0 ? ceil((x)-0.5) : floor((x)+0.5)); }
ITA_BASE_API inline float roundf(float x) { return (x<0 ? ceil((x)-0.5f) : floor((x)+0.5f)); }
#endif

//! Nächstkleineres Vielfaches von mul zu einer Zahl x zurückgeben
/**
 * Beispiele: lwrmul(10, 5) = 10
 *            lwrmul(10, 6) = 6
 *            lwrmul(-10, 5) = -10
 *            lwrmul(-10, -6) = -12 (Kleiner!)
 *            lwrmul(-10, 6) = -12
 *            lwrmul(10, -6) = 12
 */
ITA_BASE_API int lwrmul(int x, int mul);

//! Nächstgrößeres Vielfaches von mul zu einer Zahl x zurückgeben
/**
 * Beispiele: uprmul(10, 5) = 10
 *            uprmul(10, 6) = 12
 *            uprmul(-10, 5) = -10
 *            uprmul(-10, -6) = -6 (Größer!)
 *            uprmul(-10, 6) = -6
 *            uprmul(10, -6) = 12
 */
ITA_BASE_API int uprmul(int x, int mul);

//! Nächstkleineres Vielfaches von mul zu einer Zahl x zurückgeben
/**
 * Beispiele: lwrmulu(10, 5) = 10
 *            lwrmulu(10, 6) = 6
 */
ITA_BASE_API unsigned int lwrmulu(unsigned int x, unsigned int mul);

//! Nächstgrößeres Vielfaches von mul zu einer Zahl x zurückgeben
/**
 * Beispiele: uprmulu(10, 5) = 10
 *            uprmulu(10, 6) = 12
 */
ITA_BASE_API unsigned int uprmulu(unsigned int x, unsigned int mul);

//! Aufrundende Ganzzahl-Division
/**
 * Fährt die Ganzzahl-Division c=a/b durch. Falls b nicht a teil
 * (d.h. es bleibt ein Rest a mod b != 0 übrig), so wird c um
 * eins inkrementiert. Danach gibt die Funktion c zurück
 *
 * Anwendungen: Wieviele Blöcke der Länge b sind erforderlich 
 *              um ein Signal der Länge l zu beschreiben: uprdiv(l, b)
 */
ITA_BASE_API int uprdiv(int a, int b);

//! Aufrundende Ganzzahl-Division
/**
 * Variante von uprdiv für vorzeichenlose Ganzzahlen.
 */
ITA_BASE_API unsigned int uprdivu(unsigned int a, unsigned int b);

//! Einen (beliebigen) Winkel ins Interval (-180,180] abbilden
/**
 * Beispiel: correctAngle(380°) = 20°
 * \note Der exakte Wert -180 wird durch die Funktion erhalten
 *       und nicht auf +180 abgebildet (Ästethik :-))
 */
ITA_BASE_API float correctAngle180(float phi);

//! Einen (beliebigen) Winkel ins Interval [0,360) abbilden
/**
 * Beispiel: correctAngle(380°) = 20°
 */
ITA_BASE_API float correctAngle360(float phi);

//! Verhältnis in Dezibel (Energie) (im 10er-Logarithmus) umrechnen
/**
 * \param db Dezibel (Wertebereich: -Unendlich..+Unendlich)
 * \return Verhältnis (Wertebereich: 0..+Unendlich)
 *
 * \note Für die symbolische Konstante DECIBEL_MINUS_INFINITY gibt die
 *       Funktion den Wert 0 zurück.
 */
ITA_BASE_API double db10_to_ratio(double db);

//! Dezibel (Energie) (im 10er-Logarithmus) in Verhältnis umrechnen
/**
 * \param Verhältnis (Wertebereich: 0..+Unendlich)
 * \return db Dezibel (Wertebereich: -Unendlich..+Unendlich)
 *
 * \note Für den Wert 0 gibt die Funktion die die symbolische Konstante
 *       DECIBEL_MINUS_INFINITY zurück.
 *
 * \important Negative Verhältnisse sind nicht erlaubt. In diesem
 *            Falle wird eine ITAException (INVALID_PARAMETER) ausgelöst!
 */
ITA_BASE_API double ratio_to_db10(double r);


//! Verhältnis in Dezibel (Spannung) (im 10er-Logarithmus) umrechnen
/**
 * \param db Dezibel (Wertebereich: -Unendlich..+Unendlich)
 * \return Verhältnis (Wertebereich: 0..+Unendlich)
 *
 * \note Für die symbolische Konstante DECIBEL_MINUS_INFINITY gibt die
 *       Funktion den Wert 0 zurück.
 */
ITA_BASE_API double db20_to_ratio(double db);

//! Dezibel (Spannung) (im 10er-Logarithmus) in Verhältnis umrechnen
/**
 * \param Verhältnis (Wertebereich: 0..+Unendlich)
 * \return db Dezibel (Wertebereich: -Unendlich..+Unendlich)
 *
 * \note Für den Wert 0 gibt die Funktion die die symbolische Konstante
 *       DECIBEL_MINUS_INFINITY zurück.
 *
 * \important Negative Verhältnisse sind nicht erlaubt. In diesem
 *            Falle wird eine ITAException (INVALID_PARAMETER) ausgelöst!
 */
ITA_BASE_API double ratio_to_db20(double r);

/* +------------------------------------+
   |                                    |
   |   Funktionen für komplexe Zahlen   |
   |                                    |
   +------------------------------------+ */

//! Betrag einer komplexen Zahl berechnen
176
ITA_BASE_API float cabsf( float re, float im );
Jonas Stienen's avatar
Jonas Stienen committed
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429

//! Phase einer komplexen Zahl berechnen
ITA_BASE_API float canglef(float re, float im);

//! Betrag einer komplexen Zahl setzen, deren Phase aber erhalten (auch in-place)
ITA_BASE_API void csabsparg(const float& in_re, const float& in_im, const float& dest_abs, float& out_re, float& out_im);

//! Winkel einer komplexen Zahl setzen, deren Betrag aber erhalten (auch in-place)
ITA_BASE_API void csargpabs(const float& in_re, const float& in_im, const float& dest_arg, float& out_re, float& out_im);

/* +---------------------------+
   |                           |
   |   Neue Winkelfunktionen   |
   |                           |
   +---------------------------+ */

//! Winkel im Bogenmaß in Grad [°] umrechnen
ITA_BASE_API float rad2gradf(float phi);

//! Winkel im Bogenmaß in Grad [°] umrechnen
ITA_BASE_API double rad2grad(double phi);

//! Winkel in Grad [°] ins Bogenmaß umrechnen
ITA_BASE_API float grad2radf(float phi);

//! Winkel im Bogenmaß in Grad [°] umrechnen
ITA_BASE_API double grad2rad(double phi);


/* ---------- RAD ----------- */

//! Einen Winkel (rad) in das Intervall [0, 2PI) projezieren
/**
 * Projeziert einen Winkel in Rad [rad] (nicht Grad) in das Intervall [0,2*PI) und
 * setzt somit negative Winkel (gegen den UZS) in positive (im UZS) um und entfernt 
 * Mehrfachumdrehungen.
 *
 * Beispiele: 0->0, 2*PI->0, (-2*PI)->0, (2*PI+0.1)->0.1, (-0.1)->(2*PI-0.1)
 */
ITA_BASE_API float anglef_proj_0_2PI(float alpha);

//! Einen Winkel (rad) in das Intervall (-PI, PI] projezieren
/**
 * Projeziert einen Winkel in Rad [rad] (nicht Grad) in das Intervall (-PI,PI] und
 * setzt somit negative Winkel über -PI (gegen den UZS) in positive (im UZS) um sowie
 * positive Winkel über PI (im UZS) in negative (gegen UZS) und entfernt Mehrfachumdrehungen.
 * Bemerkung: -PI wird auf PI abgebildet.
 *
 * Beispiele: 0->0, 2*PI->0, (-2*PI)->0, (PI+0.1)->(-PI+0.1), (-PI-0.1)->(PI-0.1)
 */
ITA_BASE_API float anglef_proj_NPI_PI(float alpha);

//! (Gerichtete) minimale Winkeldifferenz für zwei Winkel [rad] (nicht Grad) im Intervall [0,2PI)
/**
 * Für zwei gegebene Winkel 0<=alpha, beta<2*PI im Interval [0,2*PI) berechnet
 * die Funktion den betragskleinsten Winkel gamma, um den der Winkel alpha (weiter)gedreht
 * werden muss, um den Wert des Winkels beta zu erreichen. Dies kann durch Linksdrehung
 * (im UZS, positiver Drehwinkel gamma) oder Rechtsdrehung (gegen UZS, negativer Drehwinkel gamma)
 * geschehen. Die Funktion berechnet den betragskleinsten Drehwinkel (im UZS oder gegen den UZS)
 * und gibt ihn zurück. Das Ergebnis ist also die Drehung um den geringstmöglichen Winkel um
 * alpha in beta zu überführen.
 *
 * \return Minimaler Differenzwinkel (negatives Vorzeichen = Linksdrehung, positives Vorzeichen = Rechtsdrehung)
 *
 * \note Im Gegensatz zur Funktion anglef_mindiff_abs_0_2PI gibt hier das Vorzeichen die
 *       Drehrichtung an.
 *
 * \note Die Funktion ist nicht kommutativ: f(a,b) != f(b,a)
 *
 * Beispiele: (0,PI/2) -> PI/2, (PI/5,PI/4) -> PI/20, (PI/3,PI/4) -> -PI/12, (0,2*PI-0.1) -> -0.1
 */
ITA_BASE_API float anglef_mindiff_0_2PI(float alpha, float beta);

//! Absolute minimale Winkeldifferenz für zwei Winkel [rad] im Intervall [0,2PI)
/**
 * Für zwei gegebene Winkel 0<=alpha, beta<2*PI im Interval [0,2*PI) berechnet
 * die Funktion den minimalen Differenzwinkel. Dieser Winkel ist der kleinste Winkel
 * um den einer der beiden Winkel rotiert werden muss, um den Wert des anderen Winkels
 * zu erreichen. Hierbei wird keine Drehrichtung beachtet und der Rückgabewert ist
 * stets positiv.
 *
 * \note Die Funktion ist kommutativ: f(a,b) = f(b,a)
 *
 * Beispiele: (0,PI/2) -> PI/2, (PI/2,PI/3) -> PI/6, (PI/3,PI/2) -> PI/6, (0,2*PI-0.1) -> 0.1
 */
ITA_BASE_API float anglef_mindiff_abs_0_2PI(float alpha, float beta);

/* ---------- GRAD ---------- */

//! Einen Winkel [°] in das Intervall [0°,360°) projezieren
/**
 * Projeziert einen Winkel in Grad [°] (nicht Bogenmaß) in das Intervall [0°,360°) und
 * setzt somit negative Winkel (gegen den UZS) in positive (im UZS) um und entfernt 
 * Mehrfachumdrehungen.
 *
 * Beispiele: 0°->0°, 360°->0°, -360°->0°, 361°->1°, -1°->359°
 */
ITA_BASE_API float anglef_proj_0_360_DEG(float alpha);

//! Einen Winkel (°) in das Intervall (-180°, 180°] projezieren
/**
 * Projeziert einen Winkel in Grad [°] (nicht Bogenmaß) in das Intervall (-180°,180°] und
 * setzt somit negative Winkel über -PI (gegen den UZS) in positive (im UZS) um sowie
 * positive Winkel über PI (im UZS) in negative (gegen UZS) und entfernt Mehrfachumdrehungen.
 * Bemerkung: -180° wird auf 180° abgebildet.
 *
 * Beispiele: 0°->0°, 360°->0°, -360°->0°, 181°->-179°, -181°->179°
 */
ITA_BASE_API float anglef_proj_N180_180_DEG(float alpha);

//! (Gerichtete) minimale Winkeldifferenz für zwei Winkel [°] im Intervall [0°,360°)
/**
 * Für zwei gegebene Winkel 0°<= alpha, beta < 360° im Interval [0°,360°) berechnet
 * die Funktion den betragskleinsten Winkel gamma, um den der Winkel alpha (weiter)gedreht
 * werden muss, um den Wert des Winkels beta zu erreichen. Dies kann durch Linksdrehung
 * (im UZS, positiver Drehwinkel gamma) oder Rechtsdrehung (gegen UZS, negativer Drehwinkel gamma)
 * geschehen. Die Funktion berechnet den betragskleinsten Drehwinkel (im UZS oder gegen den UZS)
 * und gibt ihn zurück. Das Ergebnis ist also die Drehung um den geringstmöglichen Winkel um
 * alpha in beta zu überführen.
 *
 * \return Minimaler Differenzwinkel (negatives Vorzeichen = Linksdrehung, positives Vorzeichen = Rechtsdrehung)
 *
 * \note Im Gegensatz zur Funktion anglef_mindiff_abs_0_360_DEG gibt hier das Vorzeichen die
 *       Drehrichtung an.
 *
 * \note Die Funktion ist nicht kommutativ: f(a,b) != f(b,a)
 *
 * Beispiele: (0°,90°) -> 90°, (45°,70°) -> 25°, (70°,45°) -> -25°, (0°,359°) -> -1°
 */
ITA_BASE_API float anglef_mindiff_0_360_DEG(float alpha, float beta);

//! Absolute minimale Winkeldifferenz für zwei Winkel [°] im Intervall [0°,360°)
/**
 * Für zwei gegebene Winkel 0°<= alpha, beta < 360° im Interval [0°,360°) berechnet
 * die Funktion den minimalen Differenzwinkel. Dieser Winkel ist der kleinste Winkel
 * um den einer der beiden Winkel rotiert werden muss, um den Wert des anderen Winkels
 * zu erreichen. Hierbei wird keine Drehrichtung beachtet und der Rückgabewert ist
 * stets positiv.
 *
 * \note Die Funktion ist kommutativ: f(a,b) = f(b,a)
 *
 * Beispiele: (0°,90°) -> 90°, (45°,70°) -> 25°, (70°,45°) -> 25°, (0°,359°) -> 1°
 */
ITA_BASE_API float anglef_mindiff_abs_0_360_DEG(float alpha, float beta);


// Orientierung in Yaw-Pitch-Roll Winkeln in View-Up-Vektor umrechnen (Alle Winkel in Bogenmaß)
ITA_BASE_API void convertYPR2VU(const double yaw, const double pitch, const double roll,
	double& vx, double& vy, double& vz, 
	double& ux, double& uy, double& uz);

// Orientierung in Yaw-Pitch-Roll Winkeln in View-Up-Vektor umrechnen (Alle Winkel in Bogenmaß)
ITA_BASE_API void convertYPR2VU(const float yaw, const float pitch, const float roll,
	float& vx, float& vy, float& vz, 
	float& ux, float& uy, float& uz);

// Orientierung in View-Up-Vektor in Yaw-Pitch-Roll Winkeln umrechnen (Alle Winkel in Bogenmaß)
ITA_BASE_API void convertVU2YPR(const double vx, const double vy, const double vz,
	const double ux, const double uy, const double uz,
	double& yaw, double& pitch, double& roll);

// Orientierung in View-Up-Vektor in Yaw-Pitch-Roll Winkeln umrechnen (Alle Winkel in Bogenmaß)
ITA_BASE_API void convertVU2YPR(const float vx, const float vy, const float vz,
	const float ux, const float uy, const float uz,
	float& yaw, float& pitch, float& roll);

// Datenklasse für Fehlerwerte
template <typename T>
class ITA_BASE_API ErrorValues {
public:
	T minAbsError;
	T avgAbsError;
	T maxAbsError;

	int indexMaxAbsError;	// Index of element where max abs error

	ErrorValues() : minAbsError(0), avgAbsError(0), maxAbsError(0), indexMaxAbsError(0) {}
};

// Fehler zwischen zwei Vektoren berechnen
template <typename Tc, typename Ta, typename Tb>
ErrorValues<Tc> computeErrorValues(const Ta* A, const Tb* B, int size)
{
	ErrorValues<Tc> v;
	if (size == 0) return v;

	Tb absErr = std::abs( (Tc) A[0] - (Tc) B[0] );
	v.minAbsError = v.avgAbsError = v.maxAbsError = absErr;

	for (int i=1; i<size; ++i) {
		absErr = std::abs( (Tc) A[i] - (Tc) B[i] );
		v.minAbsError = std::min(v.minAbsError, absErr);
		
		if (absErr > v.maxAbsError) {
			v.maxAbsError = absErr;
			v.indexMaxAbsError = i;
		}

		v.avgAbsError += absErr;
	}
	v.avgAbsError /= size;

	return v;
};

/* +-----------------------+
   |                       |
   |   Listen und Folgen   |
   |                       |
   +-----------------------+ */

// Fills a vector with numbers dest = { a+n*s < b | n in N }
template <typename T> void linspace(std::vector<T>& dest, T a, T b, T s=1) {
	dest.clear();
	if (s==0) {
		dest.push_back(a);
		return;
	}

	if (s<0) {
		if (a<b) return; // No elements contained
		dest.reserve( (a-b)/s + 1 );
		for (T x=a; x>=b; x+=s) dest.push_back(x);
	} else {
		if (a>b) return; // No elements contained
		dest.reserve( (b-a)/s + 1 );
		for (T x=a; x<=b; x+=s) dest.push_back(x);
	}
};

// Fills a vector with powers of two in the range a=2^i <= 2^j <= 2^k=b
// (Note a and b must be powers of two, otherwise an exception is thrown)
ITA_BASE_API void pow2space(std::vector<int>& dest, int a, int b);

// Calculates the factorial of an positive integer m
ITA_BASE_API int factorial(int m);

//Calculates the normalizing constant for SHRealvaluedBasefunctions
ITA_BASE_API double SHNormalizeConst(int m, int n);

//Calculates the Kronecker delta
ITA_BASE_API int SHKronecker(int m);

// Returns the index of a basefunction with degree m and order n
ITA_BASE_API int SHDegreeOrder2Linear(int m, int n);

//Calculates the realvalued Basefunctions of SH for e.g. Ambisonics
ITA_BASE_API std::vector<double> SHRealvaluedBasefunctions(double elevation, double azimuth, int maxOrder);

//Calculates the associated legendre polynomials
ITA_BASE_API std::vector<double> SHAssociatedLegendre(int N, double mu);

#endif // INCLUDE_WATCHER_ITA_NUMERIC_UTILS