1 /******************************************************************************\
2 * Copyright (c) 2001
3 *
4 * Author(s):
5 * Volker Fischer
6 *
7 * Description:
8 * c++ Mathematic Library (Matlib), standard toolbox
9 *
10 ******************************************************************************
11 *
12 * This program is free software; you can redistribute it and/or modify it under
13 * the terms of the GNU General Public License as published by the Free Software
14 * Foundation; either version 2 of the License, or (at your option) any later
15 * version.
16 *
17 * This program is distributed in the hope that it will be useful, but WITHOUT
18 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
20 * details.
21 *
22 * You should have received a copy of the GNU General Public License along with
23 * this program; if not, write to the Free Software Foundation, Inc.,
24 * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 \******************************************************************************/
27
28 #ifndef _MATLIB_STD_TOOLBOX_H_
29 #define _MATLIB_STD_TOOLBOX_H_
30
31 #include "Matlib.h"
32
33 /* fftw (Homepage: http://www.fftw.org) */
34 #include <fftw3.h>
35
36
37 /* Classes ********************************************************************/
38 class CFftPlans
39 {
40 public:
CFftPlans()41 CFftPlans() : RFFTPlForw(NULL), RFFTPlBackw(NULL), bInitialized(false) {}
CFftPlans(const int iFftSize)42 CFftPlans(const int iFftSize) : RFFTPlForw(NULL), RFFTPlBackw(NULL), bInitialized(false) {Init(iFftSize);}
43 virtual ~CFftPlans();
44
45 void Init(const int iFSi);
IsInitialized()46 inline bool IsInitialized() const {return bInitialized;}
47
48 fftw_plan RFFTPlForw;
49 fftw_plan RFFTPlBackw;
50 fftw_plan FFTPlForw;
51 fftw_plan FFTPlBackw;
52
53 double * pFftwRealIn;
54 double * pFftwRealOut;
55 fftw_complex* pFftwComplexIn;
56 fftw_complex* pFftwComplexOut;
57
58 protected:
59 void Clean();
60 bool bInitialized;
61 };
62
63
64 /* Helpfunctions **************************************************************/
Min(const CReal & rA,const CReal & rB)65 inline CReal Min(const CReal& rA, const CReal& rB)
66 {return rA < rB ? rA : rB;}
Min(const CMatlibVector<CReal> & rvA,const CMatlibVector<CReal> & rvB)67 inline CMatlibVector<CReal> Min(const CMatlibVector<CReal>& rvA, const CMatlibVector<CReal>& rvB)
68 {_VECOP(CReal, rvA.GetSize(), Min(rvA[i], rvB[i]));}
69 CReal Min(const CMatlibVector<CReal>& rvI);
70 void Min(CReal& rMinVal /* out */, int& iMinInd /* out */,
71 const CMatlibVector<CReal>& rvI /* in */);
72
Min(const CReal & r1,const CReal & r2,const CReal & r3,const CReal & r4)73 inline CReal Min(const CReal& r1, const CReal& r2, const CReal& r3, const CReal& r4)
74 {return Min(Min(Min(r1, r2), r3), r4);}
Min(const CReal & r1,const CReal & r2,const CReal & r3,const CReal & r4,const CReal & r5,const CReal & r6,const CReal & r7,const CReal & r8)75 inline CReal Min(const CReal& r1, const CReal& r2, const CReal& r3, const CReal& r4,
76 const CReal& r5, const CReal& r6, const CReal& r7, const CReal& r8)
77 {return Min(Min(Min(Min(Min(Min(Min(r1, r2), r3), r4), r5), r6), r7), r8);}
78
79
Max(const CReal & rA,const CReal & rB)80 inline CReal Max(const CReal& rA, const CReal& rB)
81 {return rA > rB ? rA : rB;}
Max(const CMatlibVector<CReal> & rvA,const CMatlibVector<CReal> & rvB)82 inline CMatlibVector<CReal> Max(const CMatlibVector<CReal>& rvA, const CMatlibVector<CReal>& rvB)
83 {_VECOP(CReal, rvA.GetSize(), Max(rvA[i], rvB[i]));}
84 CReal Max(const CMatlibVector<CReal>& rvI);
85 void Max(CReal& rMaxVal /* out */, int& iMaxInd /* out */,
86 const CMatlibVector<CReal>& rvI /* in */);
87
Max(const CReal & r1,const CReal & r2,const CReal & r3)88 inline CReal Max(const CReal& r1, const CReal& r2, const CReal& r3)
89 {return Max(Max(r1, r2), r3);}
Max(const CReal & r1,const CReal & r2,const CReal & r3,const CReal & r4,const CReal & r5,const CReal & r6,const CReal & r7)90 inline CReal Max(const CReal& r1, const CReal& r2, const CReal& r3, const CReal& r4,
91 const CReal& r5, const CReal& r6, const CReal& r7)
92 {return Max(Max(Max(Max(Max(Max(r1, r2), r3), r4), r5), r6), r7);}
93
94
Ones(const int iLen)95 inline CMatlibVector<CReal> Ones(const int iLen)
96 {_VECOP(CReal, iLen, (CReal) 1.0);}
Zeros(const int iLen)97 inline CMatlibVector<CReal> Zeros(const int iLen)
98 {_VECOP(CReal, iLen, (CReal) 0.0);}
99
100
Real(const CComplex & cI)101 inline CReal Real(const CComplex& cI) {return cI.real();}
Real(const CMatlibVector<CComplex> & cvI)102 inline CMatlibVector<CReal> Real(const CMatlibVector<CComplex>& cvI)
103 {_VECOP(CReal, cvI.GetSize(), Real(cvI[i]));}
104
Imag(const CComplex & cI)105 inline CReal Imag(const CComplex& cI) {return cI.imag();}
Imag(const CMatlibVector<CComplex> & cvI)106 inline CMatlibVector<CReal> Imag(const CMatlibVector<CComplex>& cvI)
107 {_VECOP(CReal, cvI.GetSize(), Imag(cvI[i]));}
108
Conj(const CComplex & cI)109 inline CComplex Conj(const CComplex& cI) {return conj(cI);}
Conj(const CMatlibVector<CComplex> & cvI)110 inline CMatlibVector<CComplex> Conj(const CMatlibVector<CComplex>& cvI)
111 {_VECOP(CComplex, cvI.GetSize(), Conj(cvI[i]));}
Conj(const CMatlibMatrix<CComplex> & cmI)112 inline CMatlibMatrix<CComplex> Conj(const CMatlibMatrix<CComplex>& cmI)
113 {_MATOP(CComplex, cmI.GetRowSize(), cmI.GetColSize(), Conj(cmI[i]));}
114
115
116 /* Absolute and angle (argument) functions */
Abs(const CReal & rI)117 inline CReal Abs(const CReal& rI) {return fabs(rI);}
Abs(const CMatlibVector<CReal> & fvI)118 inline CMatlibVector<CReal> Abs(const CMatlibVector<CReal>& fvI)
119 {_VECOP(CReal, fvI.GetSize(), Abs(fvI[i]));}
120
Abs(const CComplex & cI)121 inline CReal Abs(const CComplex& cI) {return abs(cI);}
Abs(const CMatlibVector<CComplex> & cvI)122 inline CMatlibVector<CReal> Abs(const CMatlibVector<CComplex>& cvI)
123 {_VECOP(CReal, cvI.GetSize(), Abs(cvI[i]));}
124
Angle(const CComplex & cI)125 inline CReal Angle(const CComplex& cI) {return arg(cI);}
Angle(const CMatlibVector<CComplex> & cvI)126 inline CMatlibVector<CReal> Angle(const CMatlibVector<CComplex>& cvI)
127 {_VECOP(CReal, cvI.GetSize(), Angle(cvI[i]));}
128
129
130 /* Trigonometric functions */
Sin(const CReal & fI)131 inline CReal Sin(const CReal& fI) {return sin(fI);}
132 template<class T> inline
Sin(const CMatlibVector<T> & vecI)133 CMatlibVector<T> Sin(const CMatlibVector<T>& vecI)
134 {_VECOP(T, vecI.GetSize(), sin(vecI[i]));}
135
Cos(const CReal & fI)136 inline CReal Cos(const CReal& fI) {return cos(fI);}
137 template<class T> inline
Cos(const CMatlibVector<T> & vecI)138 CMatlibVector<T> Cos(const CMatlibVector<T>& vecI)
139 {_VECOP(T, vecI.GetSize(), cos(vecI[i]));}
140
Tan(const CReal & fI)141 inline CReal Tan(const CReal& fI) {return tan(fI);}
142 template<class T> inline
Tan(const CMatlibVector<T> & vecI)143 CMatlibVector<T> Tan(const CMatlibVector<T>& vecI)
144 {_VECOP(T, vecI.GetSize(), tan(vecI[i]));}
145
Sinh(const CReal & fI)146 inline CReal Sinh(const CReal& fI) {return sinh(fI);}
147 template<class T> inline
Sinh(const CMatlibVector<T> & vecI)148 CMatlibVector<T> Sinh(const CMatlibVector<T>& vecI)
149 {_VECOP(T, vecI.GetSize(), sinh(vecI[i]));}
150
Cosh(const CReal & fI)151 inline CReal Cosh(const CReal& fI) {return cosh(fI);}
152 template<class T> inline
Cosh(const CMatlibVector<T> & vecI)153 CMatlibVector<T> Cosh(const CMatlibVector<T>& vecI)
154 {_VECOP(T, vecI.GetSize(), cosh(vecI[i]));}
155
Tanh(const CReal & fI)156 inline CReal Tanh(const CReal& fI) {return tanh(fI);}
157 template<class T> inline
Tanh(const CMatlibVector<T> & vecI)158 CMatlibVector<T> Tanh(const CMatlibVector<T>& vecI)
159 {_VECOP(T, vecI.GetSize(), tanh(vecI[i]));}
160
161
162 /* Square root */
Sqrt(const CReal & fI)163 inline CReal Sqrt(const CReal& fI) {return sqrt(fI);}
164 template<class T> inline
Sqrt(const CMatlibVector<T> & vecI)165 CMatlibVector<T> Sqrt(const CMatlibVector<T>& vecI)
166 {_VECOP(T, vecI.GetSize(), sqrt(vecI[i]));}
167
168
169 /* Exponential function */
Exp(const CReal & fI)170 inline CReal Exp(const CReal& fI) {return exp(fI);}
171 template<class T> inline
Exp(const CMatlibVector<T> & vecI)172 CMatlibVector<T> Exp(const CMatlibVector<T>& vecI)
173 {_VECOP(T, vecI.GetSize(), exp(vecI[i]));}
174
175
176 /* Logarithm */
Log(const CReal & fI)177 inline CReal Log(const CReal& fI) {return log(fI);}
178 template<class T> inline
Log(const CMatlibVector<T> & vecI)179 CMatlibVector<T> Log(const CMatlibVector<T>& vecI)
180 {_VECOP(T, vecI.GetSize(), log(vecI[i]));}
181
Log10(const CReal & fI)182 inline CReal Log10(const CReal& fI) {return log10(fI);}
183 template<class T> inline
Log10(const CMatlibVector<T> & vecI)184 CMatlibVector<T> Log10(const CMatlibVector<T>& vecI)
185 {_VECOP(T, vecI.GetSize(), log10(vecI[i]));}
186
187
188 /* Mean, variance and standard deviation */
Mean(const CMatlibVector<T> & vecI)189 template<class T> inline T Mean(const CMatlibVector<T>& vecI)
190 {return Sum(vecI) / vecI.GetSize();}
Std(CMatlibVector<T> & vecI)191 template<class T> inline T Std(CMatlibVector<T>& vecI)
192 {return Sqrt(Var(vecI));}
193 template<class T> T Var(const CMatlibVector<T>& vecI);
194
195
196 /* Rounding functions */
Fix(const CReal & fI)197 inline CReal Fix(const CReal& fI) {return (int) fI;}
Fix(const CMatlibVector<CReal> & fvI)198 inline CMatlibVector<CReal> Fix(const CMatlibVector<CReal>& fvI)
199 {_VECOP(CReal, fvI.GetSize(), Fix(fvI[i]));}
200
Floor(const CReal & fI)201 inline CReal Floor(const CReal& fI) {return floor(fI);}
Floor(const CMatlibVector<CReal> & fvI)202 inline CMatlibVector<CReal> Floor(const CMatlibVector<CReal>& fvI)
203 {_VECOP(CReal, fvI.GetSize(), Floor(fvI[i]));}
204
Ceil(const CReal & fI)205 inline CReal Ceil(const CReal& fI) {return ceil(fI);}
Ceil(const CMatlibVector<CReal> & fvI)206 inline CMatlibVector<CReal> Ceil(const CMatlibVector<CReal>& fvI)
207 {_VECOP(CReal, fvI.GetSize(), Ceil(fvI[i]));}
208
Round(const CReal & fI)209 inline CReal Round(const CReal& fI)
210 {return Floor(fI + (CReal) 0.5);}
Round(const CMatlibVector<CReal> & fvI)211 inline CMatlibVector<CReal> Round(const CMatlibVector<CReal>& fvI)
212 {_VECOP(CReal, fvI.GetSize(), Round(fvI[i]));}
213
Sign(const CReal & rI)214 inline CReal Sign(const CReal& rI)
215 {return rI == 0 ? 0 : rI > 0 ? 1 : -1;}
216
Mod(const int ix,const int iy)217 inline int Mod(const int ix, const int iy)
218 {return ix < 0 ? (ix % iy + iy) % iy : ix % iy;}
219
220 template<class T> T Sum(const CMatlibVector<T>& vecI);
221
222 CMatlibVector<CReal> Sort(const CMatlibVector<CReal>& rvI);
223
224
225 /* Matrix inverse */
226 CMatlibMatrix<CComplex> Inv(const CMatlibMatrix<CComplex>& matrI);
227
228 /* Identity matrix */
229 CMatlibMatrix<CReal> Eye(const int iLen);
230
231 CMatlibMatrix<CComplex> Diag(const CMatlibVector<CComplex>& cvI);
232
233 CReal Trace(const CMatlibMatrix<CReal>& rmI);
234
235 CMatlibMatrix<CComplex> Toeplitz(const CMatlibVector<CComplex>& cvI);
236
237 /* Matrix transpose */
238 CMatlibMatrix<CComplex> Transp(const CMatlibMatrix<CComplex>& cmI);
239 inline
TranspH(const CMatlibMatrix<CComplex> & cmI)240 CMatlibMatrix<CComplex> TranspH(const CMatlibMatrix<CComplex>& cmI)
241 {return Conj(Transp(cmI));} /* With conjugate complex */
242
243 /* Fourier transformations (also included: real FFT) */
244 CMatlibVector<CComplex> Fft(const CMatlibVector<CComplex>& cvI, const CFftPlans& FftPlans = CFftPlans());
245 CMatlibVector<CComplex> Ifft(const CMatlibVector<CComplex>& cvI, const CFftPlans& FftPlans = CFftPlans());
246 CMatlibVector<CComplex> rfft(const CMatlibVector<CReal>& fvI, const CFftPlans& FftPlans = CFftPlans());
247 CMatlibVector<CReal> rifft(const CMatlibVector<CComplex>& cvI, const CFftPlans& FftPlans = CFftPlans());
248
249 CMatlibVector<CReal> FftFilt(const CMatlibVector<CComplex>& rvH,
250 const CMatlibVector<CReal>& rvI,
251 CMatlibVector<CReal>& rvZ,
252 const CFftPlans& FftPlans = CFftPlans());
253
254
255 /* Numerical integration */
256 typedef CComplex(MATLIB_CALLBACK_QAUD)(CReal rX); /* Callback function definition */
257 CComplex Quad(MATLIB_CALLBACK_QAUD f, const CReal a,
258 const CReal b, const CReal errorBound = 1.e-6);
259
260
261 /* Implementation **************************************************************
262 (the implementation of template classes must be in the header file!) */
263 template<class T> inline
Sum(const CMatlibVector<T> & vecI)264 T Sum(const CMatlibVector<T>& vecI)
265 {
266 const int iSize = vecI.GetSize();
267 T SumRet = 0;
268 for (int i = 0; i < iSize; i++)
269 SumRet += vecI[i];
270
271 return SumRet;
272 }
273
274 template<class T> inline
Var(const CMatlibVector<T> & vecI)275 T Var(const CMatlibVector<T>& vecI)
276 {
277 const int iSize = vecI.GetSize();
278
279 /* First calculate mean */
280 T tMean = Mean(vecI);
281
282 /* Now variance (sum formula) */
283 T tRet = 0;
284 for (int i = 0; i < iSize; i++)
285 tRet += (vecI[i] - tMean) * (vecI[i] - tMean);
286
287 return tRet / (iSize - 1); /* Normalizing */
288 }
289
290
291 #endif /* _MATLIB_STD_TOOLBOX_H_ */
292