IT++ Logo
specmat.cpp
Go to the documentation of this file.
1 
30 #include <itpp/base/specmat.h>
32 #include <itpp/base/math/log_exp.h>
33 #include <itpp/base/itcompat.h>
34 #include <itpp/base/matfunc.h>
35 
36 
37 namespace itpp
38 {
39 
40 ivec find(const bvec &invector)
41 {
42  it_assert(invector.size() > 0, "find(): vector cannot be empty");
43  ivec temp(invector.size());
44  int pos = 0;
45  for (int i = 0;i < invector.size();i++) {
46  if (invector(i) == bin(1)) {
47  temp(pos) = i;
48  pos++;
49  }
50  }
51  temp.set_size(pos, true);
52  return temp;
53 }
54 
56 
57 #define CREATE_SET_FUNS(typef,typem,name,value) \
58  typef name(int size) \
59  { \
60  typef t(size); \
61  t = value; \
62  return t; \
63  } \
64  \
65  typem name(int rows, int cols) \
66  { \
67  typem t(rows, cols); \
68  t = value; \
69  return t; \
70  }
71 
72 #define CREATE_EYE_FUN(type,name,zero,one) \
73  type name(int size) { \
74  type t(size,size); \
75  t = zero; \
76  for (int i=0; i<size; i++) \
77  t(i,i) = one; \
78  return t; \
79  }
80 
81 CREATE_SET_FUNS(vec, mat, ones, 1.0)
82 CREATE_SET_FUNS(bvec, bmat, ones_b, bin(1))
83 CREATE_SET_FUNS(ivec, imat, ones_i, 1)
84 CREATE_SET_FUNS(cvec, cmat, ones_c, std::complex<double>(1.0))
85 
86 CREATE_SET_FUNS(vec, mat, zeros, 0.0)
87 CREATE_SET_FUNS(bvec, bmat, zeros_b, bin(0))
88 CREATE_SET_FUNS(ivec, imat, zeros_i, 0)
89 CREATE_SET_FUNS(cvec, cmat, zeros_c, std::complex<double>(0.0))
90 
91 CREATE_EYE_FUN(mat, eye, 0.0, 1.0)
92 CREATE_EYE_FUN(bmat, eye_b, bin(0), bin(1))
93 CREATE_EYE_FUN(imat, eye_i, 0, 1)
94 CREATE_EYE_FUN(cmat, eye_c, std::complex<double>(0.0), std::complex<double>(1.0))
95 
97 
98 vec impulse(int size)
99 {
100  vec t(size);
101  t.clear();
102  t[0] = 1.0;
103  return t;
104 }
105 
106 vec linspace(double from, double to, int points)
107 {
108  if (points < 2) {
109  // This is the "Matlab definition" of linspace
110  vec output(1);
111  output(0) = to;
112  return output;
113  }
114  else {
115  vec output(points);
116  double step = (to - from) / double(points - 1);
117  int i;
118  for (i = 0; i < (points-1); i++)
119  output(i) = from + i * step;
120  output(i) = to;
121  return output;
122  }
123 }
124 
125 vec zigzag_space(double t0, double t1, int K)
126 {
127  it_assert(K > 0, "zigzag_space:() K must be positive");
128  ivec N = "0 1";
129 
130  int n = 2;
131  for (int k = 0; k < K; k++) {
132  ivec Nn = 2 * N;
133  for (int i = 1; i < length(Nn); i += 2) {
134  Nn = concat(Nn, i);
135  n++;
136  }
137  N = Nn;
138  }
139 
140  vec T0 = linspace(t0, t1, n);
141  vec Tt = zeros(n);
142  for (int i = 0; i < n; i++) {
143  Tt(i) = T0(N(i));
144  }
145  return Tt;
146 }
147 
148 // Construct a Hadamard-imat of size "size"
149 imat hadamard(int size)
150 {
151  it_assert(size > 0, "hadamard(): size is not a power of 2");
152  int logsize = ceil_i(::log2(static_cast<double>(size)));
153  it_assert(pow2i(logsize) == size, "hadamard(): size is not a power of 2");
154 
155  imat H(size, size);
156  H(0, 0) = 1;
157 
158  for (int i = 0; i < logsize; ++i) {
159  int pow2 = 1 << i;
160  for (int k = 0; k < pow2; ++k) {
161  for (int l = 0; l < pow2; ++l) {
162  H(k, l) = H(k, l);
163  H(k + pow2, l) = H(k, l);
164  H(k, l + pow2) = H(k, l);
165  H(k + pow2, l + pow2) = (-1) * H(k, l);
166  }
167  }
168  }
169  return H;
170 }
171 
172 imat jacobsthal(int p)
173 {
174  int quadratic_residue;
175  imat out(p, p);
176  int i, j;
177 
178  out = -1; // start with all elements equal to "-1"
179 
180  // Generate a complete list of quadratic residues
181  for (i = 0; i < (p - 1) / 2; i++) {
182  quadratic_residue = ((i + 1) * (i + 1)) % p;
183  // set this element in all rows (col-row) = quadratic_residue
184  for (j = 0; j < p; j++) {
185  out(j, (j + quadratic_residue) % p) = 1;
186  }
187  }
188 
189  // set diagonal elements to zero
190  for (i = 0; i < p; i++) {
191  out(i, i) = 0;
192  }
193  return out;
194 }
195 
196 imat conference(int n)
197 {
198  it_assert_debug(n % 4 == 2, "conference(int n); wrong size");
199  int pm = n - 1; // p must be odd prime, not checked
200  imat out(n, n);
201 
202  out.set_submatrix(1, 1, jacobsthal(pm));
203  out.set_submatrix(0, 0, 1, n - 1, 1);
204  out.set_submatrix(1, n - 1, 0, 0, 1);
205  out(0, 0) = 0;
206 
207  return out;
208 }
209 
210 const cmat toeplitz(const cvec &c)
211 {
212  int s = c.size();
213  cmat output(s, s);
214  cvec c_conj = conj(c);
215  for (int i = 1; i < s; ++i) {
216  for (int j = 0; j < s - i; ++j) {
217  output(i + j, j) = c_conj(i);
218  }
219  }
220  // start from j = 0 here, because the main diagonal is not conjugated
221  for (int j = 0; j < s; ++j) {
222  for (int i = 0; i < s - j; ++i) {
223  output(i, i + j) = c(j);
224  }
225  }
226  return output;
227 }
228 
229 mat rotation_matrix(int dim, int plane1, int plane2, double angle)
230 {
231  mat m;
232  double c = std::cos(angle), s = std::sin(angle);
233 
234  it_assert(plane1 >= 0 && plane2 >= 0 &&
235  plane1 < dim && plane2 < dim && plane1 != plane2,
236  "Invalid arguments to rotation_matrix()");
237 
238  m.set_size(dim, dim, false);
239  m = 0.0;
240  for (int i = 0; i < dim; i++)
241  m(i, i) = 1.0;
242 
243  m(plane1, plane1) = c;
244  m(plane1, plane2) = -s;
245  m(plane2, plane1) = s;
246  m(plane2, plane2) = c;
247 
248  return m;
249 }
250 
251 void house(const vec &x, vec &v, double &beta)
252 {
253  double sigma, mu;
254  int n = x.size();
255 
256  v = x;
257  if (n == 1) {
258  v(0) = 1.0;
259  beta = 0.0;
260  return;
261  }
262  sigma = sum(sqr(x(1, n - 1)));
263  v(0) = 1.0;
264  if (sigma == 0.0)
265  beta = 0.0;
266  else {
267  mu = std::sqrt(sqr(x(0)) + sigma);
268  if (x(0) <= 0.0)
269  v(0) = x(0) - mu;
270  else
271  v(0) = -sigma / (x(0) + mu);
272  beta = 2 * sqr(v(0)) / (sigma + sqr(v(0)));
273  v /= v(0);
274  }
275 }
276 
277 void givens(double a, double b, double &c, double &s)
278 {
279  double t;
280 
281  if (b == 0) {
282  c = 1.0;
283  s = 0.0;
284  }
285  else {
286  if (fabs(b) > fabs(a)) {
287  t = -a / b;
288  s = -1.0 / std::sqrt(1 + t * t);
289  c = s * t;
290  }
291  else {
292  t = -b / a;
293  c = 1.0 / std::sqrt(1 + t * t);
294  s = c * t;
295  }
296  }
297 }
298 
299 void givens(double a, double b, mat &m)
300 {
301  double t, c, s;
302 
303  m.set_size(2, 2);
304 
305  if (b == 0) {
306  m(0, 0) = 1.0;
307  m(1, 1) = 1.0;
308  m(1, 0) = 0.0;
309  m(0, 1) = 0.0;
310  }
311  else {
312  if (fabs(b) > fabs(a)) {
313  t = -a / b;
314  s = -1.0 / std::sqrt(1 + t * t);
315  c = s * t;
316  }
317  else {
318  t = -b / a;
319  c = 1.0 / std::sqrt(1 + t * t);
320  s = c * t;
321  }
322  m(0, 0) = c;
323  m(1, 1) = c;
324  m(0, 1) = s;
325  m(1, 0) = -s;
326  }
327 }
328 
329 mat givens(double a, double b)
330 {
331  mat m(2, 2);
332  givens(a, b, m);
333  return m;
334 }
335 
336 void givens_t(double a, double b, mat &m)
337 {
338  double t, c, s;
339 
340  m.set_size(2, 2);
341 
342  if (b == 0) {
343  m(0, 0) = 1.0;
344  m(1, 1) = 1.0;
345  m(1, 0) = 0.0;
346  m(0, 1) = 0.0;
347  }
348  else {
349  if (fabs(b) > fabs(a)) {
350  t = -a / b;
351  s = -1.0 / std::sqrt(1 + t * t);
352  c = s * t;
353  }
354  else {
355  t = -b / a;
356  c = 1.0 / std::sqrt(1 + t * t);
357  s = c * t;
358  }
359  m(0, 0) = c;
360  m(1, 1) = c;
361  m(0, 1) = -s;
362  m(1, 0) = s;
363  }
364 }
365 
366 mat givens_t(double a, double b)
367 {
368  mat m(2, 2);
369  givens_t(a, b, m);
370  return m;
371 }
372 
374 template void eye(int, mat &);
376 template void eye(int, bmat &);
378 template void eye(int, imat &);
380 template void eye(int, cmat &);
381 
383 template vec linspace_fixed_step(double, double, double);
385 template ivec linspace_fixed_step(int, int, int);
387 template svec linspace_fixed_step(short int, short int, short int);
388 
389 } // namespace itpp
Various functions on vectors and matrices - header file.
ITPP_EXPORT ivec zeros_i(int size)
A Int vector of zeros.
ITPP_EXPORT bvec ones_b(int size)
A Binary vector of ones.
ITPP_EXPORT ivec ones_i(int size)
A Int vector of ones.
ITPP_EXPORT cmat eye_c(int size)
A Double Complex (size,size) unit matrix.
int ceil_i(double x)
The nearest larger integer.
Definition: converters.h:339
vec impulse(int size)
Impulse vector.
Definition: specmat.cpp:98
void givens_t(double a, double b, mat &m)
Calculate the transposed Givens rotation matrix.
Definition: specmat.cpp:336
ivec find(const bvec &invector)
Return a integer vector with indicies where bvec == 1.
Definition: specmat.cpp:40
cvec conj(const cvec &x)
Conjugate of complex value.
Definition: elem_math.h:226
ITPP_EXPORT bvec zeros_b(int size)
A Binary vector of zeros.
STL namespace.
template void eye(int, mat &)
Template instantiation of eye.
T sum(const Vec< T > &v)
Sum of all elements in the vector.
Definition: matfunc.h:59
#define it_assert(t, s)
Abort if t is not true.
Definition: itassert.h:94
Logarithmic and exponenential functions - header file.
int length(const Vec< T > &v)
Length of vector.
Definition: matfunc.h:51
#define it_assert_debug(t, s)
Abort if t is not true and NDEBUG is not defined.
Definition: itassert.h:107
ITPP_EXPORT cvec zeros_c(int size)
A Double Complex vector of zeros.
const cmat toeplitz(const cvec &c)
Generate symmetric Toeplitz matrix from vector c (complex valued)
Definition: specmat.cpp:210
vec sin(const vec &x)
Sine function.
Definition: trig_hyp.h:54
ITPP_EXPORT vec zeros(int size)
A Double vector of zeros.
vec angle(const cvec &x)
Angle.
Definition: elem_math.h:218
vec log2(const vec &x)
log-2 of the elements
Definition: log_exp.cpp:36
int pow2i(int x)
Calculate two to the power of x (2^x); x is integer.
Definition: log_exp.h:53
double pow2(double x)
Calculate two to the power of x (2^x)
Definition: log_exp.h:55
int size(const Vec< T > &v)
Length of vector.
Definition: matfunc.h:55
itpp namespace
Definition: itmex.h:36
imat conference(int n)
Conference matrix.
Definition: specmat.cpp:196
imat jacobsthal(int p)
Jacobsthal matrix.
Definition: specmat.cpp:172
ITPP_EXPORT vec ones(int size)
A float vector of ones.
mat rotation_matrix(int dim, int plane1, int plane2, double angle)
Create a rotation matrix that rotates the given plane angle radians. Note that the order of the plane...
Definition: specmat.cpp:229
void house(const vec &x, vec &v, double &beta)
Calcualte the Householder vector.
Definition: specmat.cpp:251
IT++ compatibility types and functions.
vec sqr(const cvec &data)
Absolute square of elements.
Definition: elem_math.cpp:36
Definitions of special vectors and matrices.
template vec linspace_fixed_step(double, double, double)
Template instantiation of linspace_fixed_step.
vec linspace(double from, double to, int points)
linspace (works in the same way as the MATLAB version)
Definition: specmat.cpp:106
Binary arithmetic (boolean) class.
Definition: binary.h:56
vec sqrt(const vec &x)
Square root of the elements.
Definition: elem_math.h:123
T to(double x)
Convert double to T.
vec zigzag_space(double t0, double t1, int K)
Zig-zag space function (variation on linspace)
Definition: specmat.cpp:125
ITPP_EXPORT bmat eye_b(int size)
A Binary (size,size) unit matrix.
vec cos(const vec &x)
Cosine function.
Definition: trig_hyp.h:58
ITPP_EXPORT cvec ones_c(int size)
A float Complex vector of ones.
ITPP_EXPORT imat eye_i(int size)
A Int (size,size) unit matrix.
Elementary mathematical functions - header file.
void givens(double a, double b, double &c, double &s)
Calculate the Givens rotation values.
Definition: specmat.cpp:277
imat hadamard(int size)
Hadamard matrix.
Definition: specmat.cpp:149
Mat< bin > bmat
bin matrix
Definition: mat.h:508
const Array< T > concat(const Array< T > &a, const T &e)
Append element e to the end of the Array a.
Definition: array.h:486
SourceForge Logo

Generated on Sun Apr 10 2022 12:00:00 for IT++ by Doxygen 1.8.14