xue@11
|
1 /*
|
xue@11
|
2 Harmonic sinusoidal modelling and tools
|
xue@11
|
3
|
xue@11
|
4 C++ code package for harmonic sinusoidal modelling and relevant signal processing.
|
xue@11
|
5 Centre for Digital Music, Queen Mary, University of London.
|
xue@11
|
6 This file copyright 2011 Wen Xue.
|
xue@11
|
7
|
xue@11
|
8 This program is free software; you can redistribute it and/or
|
xue@11
|
9 modify it under the terms of the GNU General Public License as
|
xue@11
|
10 published by the Free Software Foundation; either version 2 of the
|
xue@11
|
11 License, or (at your option) any later version.
|
xue@11
|
12 */
|
xue@1
|
13 #ifndef XCOMPLEX
|
xue@1
|
14 #define XCOMPLEX
|
xue@1
|
15
|
xue@1
|
16
|
Chris@5
|
17 #include <math.h>
|
Chris@5
|
18
|
Chris@5
|
19 /**
|
Chris@5
|
20 \file xcomplex.h - Xue's complex number class.
|
xue@1
|
21
|
xue@1
|
22 This classes is modelled after standard c++ complex class, with a few additional methods. Unused
|
xue@1
|
23 standard member and non-member functions/operators may not have been implemented.
|
xue@1
|
24 */
|
xue@1
|
25
|
xue@1
|
26 template <class T> class cmplx
|
xue@1
|
27 {
|
xue@1
|
28 public:
|
xue@1
|
29 //standard members
|
xue@1
|
30 T x;
|
xue@1
|
31 T y;
|
xue@1
|
32 cmplx(){}
|
xue@1
|
33 cmplx(const T& c, const T& d) {x=c; y=d;}
|
xue@1
|
34 cmplx(const T& c) {x=c; y=0;}
|
xue@1
|
35 template<class X> cmplx(const cmplx<X>& c){x=c.x, y=c.y;}
|
xue@1
|
36
|
xue@1
|
37 T real(){return x;}
|
xue@1
|
38 T imag(){return y;}
|
xue@1
|
39
|
xue@1
|
40 cmplx<T>& operator=(const T& c){x=c; y=0; return *this;}
|
xue@10
|
41 template<class X> cmplx<T>& operator+=(const X& c){x+=c; return *this;}
|
xue@10
|
42 template<class X> cmplx<T>& operator-=(const X& c){x-=c; return *this;}
|
xue@1
|
43 template<class X> cmplx<T>& operator*=(const X& c){x*=c; y*=c; return *this;}
|
xue@1
|
44 template<class X> cmplx<T>& operator/=(const X& c){x/=c; y/=c; return *this;}
|
xue@1
|
45 template<class X> cmplx<T>& operator=(const cmplx<X>& c){x=c.x; y=c.y; return* this;}
|
xue@1
|
46
|
xue@1
|
47 template<class X> cmplx<T>& operator+=(const cmplx<X>& c){x+=c.x; y+=c.y; return *this;}
|
xue@1
|
48 template<class X> cmplx<T>& operator-=(const cmplx<X>& c){x-=c.x; y-=c.y; return *this;}
|
xue@1
|
49 template<class X> cmplx<T>& operator*=(const cmplx<X>& c){T tmpx=x*c.x-y*c.y; y=x*c.y+y*c.x; x=tmpx; return *this;}
|
xue@1
|
50 template<class X> cmplx<T>& operator/=(const cmplx<X>& c){if (c.y==0){x/=c.x; y/=c.x;} else if (c.x==0){T tmpx=y/c.y; y=-x/c.y; x=tmpx;} else {T xm=c.x*c.x+c.y*c.y; T tmpx=(c.x*x+c.y*y)/xm; y=(y*c.x-x*c.y)/xm; x=tmpx;} return *this;}
|
xue@1
|
51
|
xue@1
|
52 //non-standard members
|
xue@1
|
53 //operator~: square of absolute value
|
xue@1
|
54 T operator~() {return x*x+y*y;}
|
xue@1
|
55 //operator*: complex conjugate
|
xue@1
|
56 cmplx<T> operator*(){cmplx<T> result; result.x=x; result.y=-y; return result;}
|
xue@1
|
57 //operator^: multiplicaiton with the complex conjugate of the argument
|
xue@1
|
58 cmplx<T> operator^(const cmplx &b){cmplx<T> result; result.x=x*b.x+y*b.y; result.y=y*b.x-x*b.y; return result;}
|
xue@1
|
59 //cinv: complex reciprocal
|
xue@1
|
60 cmplx<T> cinv(){cmplx<T> result; if (y==0) result.x=1/x, result.y=0; else if (x==0) result.y=-1/y, result.x=0; else{T xm=x*x+y*y; result.x=x/xm; result.y=-y/xm;} return result;}
|
xue@1
|
61 //rotate: rotate by an angle
|
xue@1
|
62 cmplx<T>& rotate(const T ph) {double s=sin(ph), c=cos(ph), tmpx; tmpx=x*c-y*s; y=x*s+y*c; x=tmpx; return *this;}
|
xue@1
|
63 }; //*/
|
xue@1
|
64 //standard Non-member Operators
|
xue@1
|
65 template<class Ta, class Tb> cmplx<Ta> operator+(const cmplx<Ta>& a, const cmplx<Tb>& b){cmplx<Ta> result=a; result+=b; return result;}
|
xue@1
|
66 template<class Ta, class Tb> cmplx<Ta> operator+(const cmplx<Ta>& a, Tb& b){cmplx<Ta> result=a; result+=b; return result;}
|
xue@1
|
67 template<class T> cmplx<T> operator+(T a, const cmplx<T>& b){cmplx<T> result=a; result+=b; return result;}
|
xue@1
|
68 template<class Ta, class Tb> cmplx<Ta> operator-(const cmplx<Ta>& a, const cmplx<Tb>& b){cmplx<Ta> result=a; result-=b; return result;}
|
xue@1
|
69 template<class Ta, class Tb> cmplx<Ta> operator-(const cmplx<Ta>& a, Tb& b){cmplx<Ta> result=a; result-=b; return result;}
|
xue@1
|
70 template<class T> cmplx<T> operator-(T a, const cmplx<T>& b){cmplx<T> result=a; result-=b; return result;}
|
xue@1
|
71 template<class Ta, class Tb> cmplx<Ta> operator*(const cmplx<Ta>& a, const cmplx<Tb>& b){cmplx<Ta> result=a; result*=b; return result;}
|
xue@1
|
72 template<class Ta, class Tb> cmplx<Ta> operator*(const cmplx<Ta>& a, Tb& b){cmplx<Ta> result=a; result*=b; return result;}
|
xue@1
|
73 template<class T> cmplx<T> operator*(T& a, const cmplx<T>& b){cmplx<T> result=b; result*=a; return result;}
|
xue@1
|
74 template<class Ta, class Tb> cmplx<Ta> operator/(const cmplx<Ta>& a, const cmplx<Tb>& b){cmplx<Ta> result=a; result/=b; return result;}
|
xue@1
|
75 template<class Ta, class Tb> cmplx<Ta> operator/(const cmplx<Ta>& a, Tb& b){cmplx<Ta> result=a; result/=b; return result;}
|
xue@1
|
76 template<class T> cmplx<T> operator/(T a, const cmplx<T>& b){cmplx<T> result=a; result/=b; return result;}
|
xue@1
|
77 template<class T> cmplx<T> operator+(const cmplx<T>& a){return a;}
|
xue@1
|
78 template<class T> cmplx<T> operator-(const cmplx<T>& a){cmplx<T> result; result.x=-a.x; result.y=-a.y; return result;}
|
xue@1
|
79 template<class Ta, class Tb> bool operator==(const cmplx<Ta>& a, const cmplx<Tb>& b){return (a.x==b.x && a.y==b.y);}
|
xue@1
|
80 template<class Ta, class Tb> bool operator==(const cmplx<Ta>& a, Tb b){return (a.x==b && a.y==0);}
|
xue@1
|
81 template<class Ta, class Tb> bool operator==(Ta a, const cmplx<Tb>& b){return (a==b.x && 0==b.y);}
|
xue@1
|
82 template<class Ta, class Tb> bool operator!=(const cmplx<Ta>& a, const cmplx<Tb>& b){return (a.x!=b.x || a.y!=b.y);}
|
xue@1
|
83 template<class Ta, class Tb> bool operator!=(const cmplx<Ta>& a, Tb b){return (a.x!=b || a.y!=0);}
|
xue@1
|
84 template<class Ta, class Tb> bool operator!=(Ta a, const cmplx<Tb>& b){return (a!=b.x || 0!=b.y);}
|
xue@1
|
85 /*
|
xue@1
|
86 template <class T, class charT, class traits> basic_istream<charT, traits>& operator>>(istream&, complex<T>&);
|
xue@1
|
87 template <class T, class charT, class traits> basic_ostream<charT, traits>& operator<<(ostream&, const complex<T>&);
|
xue@1
|
88 */
|
xue@1
|
89 //Values
|
xue@1
|
90 template<class T> T real(const cmplx<T>& a){return a.x;}
|
xue@1
|
91 template<class T> T imag(const cmplx<T>& a){return a.y;}
|
xue@1
|
92 template<class T> T abs(const cmplx<T>& a){return sqrt(a.x*a.x+a.y*a.y);}
|
xue@1
|
93 template<class T> T fabs(const cmplx<T>& a){return sqrt(a.x*a.x+a.y*a.y);}
|
xue@1
|
94 template<class T> T arg(const cmplx<T>& a){return (a.x==0 && a.y==0)?0:atan2(a.y, a.x);}
|
xue@1
|
95 template<class T> T norm(const cmplx<T>& a){return a.x*a.x+a.y*a.y;}
|
xue@1
|
96 template<class T> cmplx<T> conj(const cmplx<T>& a){cmplx<T> result; result.x=a.x; result.y=-a.y; return result;}
|
xue@1
|
97 template<class T> cmplx<T> polar(const T& r, const T& theta){cmplx<T> result; result.x=r*cos(theta); result.y=r*sin(theta); return result;}
|
xue@1
|
98 //Transcendentals
|
xue@1
|
99 /*
|
xue@1
|
100 template<class T> cmplx<T> cos (const cmplx<T>&);
|
xue@1
|
101 template<class T> cmplx<T> cosh (const cmplx<T>&);
|
xue@1
|
102 */
|
xue@1
|
103 template<class T> cmplx<T> exp(const cmplx<T>& a){return polar(exp(a.x), a.y);}
|
xue@1
|
104 template<class T> cmplx<T> log(const cmplx<T>& a){cmplx<T> result; result.x=0.5*log(norm(a)); result.y=arg(a); return result;}
|
xue@1
|
105 /*
|
xue@1
|
106 template<class T> cmplx<T> log10 (const cmplx<T>&);
|
xue@1
|
107 template<class T> cmplx<T> pow (const cmplx<T>&, int);
|
xue@1
|
108 */
|
xue@1
|
109 template<class T> cmplx<T> pow(const cmplx<T>& a, T& e){cmplx<T> result; T rad=abs(a); T angle=arg(a); return polar(rad, angle*e);}
|
xue@1
|
110 /*
|
xue@1
|
111 template<class T> cmplx<T> pow (const cmplx<T>&, const cmplx<T>&);
|
xue@1
|
112 template<class T> cmplx<T> pow (const T&, const cmplx<T>&);
|
xue@1
|
113 template<class T> cmplx<T> sin (const cmplx<T>&);
|
xue@1
|
114 template<class T> cmplx<T> sinh (const cmplx<T>&);
|
xue@1
|
115 */
|
xue@1
|
116 template<class T> cmplx<T> sqrt(const cmplx<T>& a){cmplx<T> result; if (a.y==0) {if (a.x>=0){result.x=sqrt(a.x); result.y=0;} else{result.y=sqrt(-a.x); result.x=0;}} else {result.y=sqrt((sqrt(a.x*a.x+a.y*a.y)-a.x)*0.5); result.x=0.5*a.y/result.y;} return result;}
|
xue@1
|
117 /*
|
xue@1
|
118 template<class T> cmplx<T> tan (const cmplx<T>&);
|
xue@1
|
119 template<class T> cmplx<T> tanh (const cmplx<T>&);
|
xue@1
|
120 */
|
xue@1
|
121
|
xue@1
|
122 //non-standard non-member functions
|
xue@1
|
123 //template operator^: multiplying one complex number with the complex conjugate of another
|
xue@1
|
124 template<class T> cmplx<T> operator^(const cmplx<T>& a, const cmplx<T>& b){cmplx<T> result=a; result^=b; return result;}
|
xue@1
|
125
|
xue@1
|
126 typedef cmplx<double> cdouble;
|
xue@1
|
127 typedef cmplx<float> cfloat;
|
xue@1
|
128
|
xue@1
|
129 #endif
|