specmat.cpp
Go to the documentation of this file.
#include <itpp/base/specmat.h>
#include <itpp/base/math/elem_math.h>
#include <itpp/base/math/log_exp.h>
#include <itpp/base/itcompat.h>
#include <itpp/base/matfunc.h>
namespace itpp
{
ivec find(const bvec &invector)
{
it_assert(invector.size() > 0, "find(): vector cannot be empty");
ivec temp(invector.size());
int pos = 0;
for (int i = 0;i < invector.size();i++) {
if (invector(i) == bin(1)) {
temp(pos) = i;
pos++;
}
}
temp.set_size(pos, true);
return temp;
}
#define CREATE_SET_FUNS(typef,typem,name,value) \
typef name(int size) \
{ \
typef t(size); \
t = value; \
return t; \
} \
\
typem name(int rows, int cols) \
{ \
typem t(rows, cols); \
t = value; \
return t; \
}
#define CREATE_EYE_FUN(type,name,zero,one) \
type name(int size) { \
type t(size,size); \
t = zero; \
for (int i=0; i<size; i++) \
t(i,i) = one; \
return t; \
}
CREATE_SET_FUNS(vec, mat, ones, 1.0)
CREATE_SET_FUNS(bvec, bmat, ones_b, bin(1))
CREATE_SET_FUNS(ivec, imat, ones_i, 1)
CREATE_SET_FUNS(cvec, cmat, ones_c, std::complex<double>(1.0))
CREATE_SET_FUNS(vec, mat, zeros, 0.0)
CREATE_SET_FUNS(bvec, bmat, zeros_b, bin(0))
CREATE_SET_FUNS(ivec, imat, zeros_i, 0)
CREATE_SET_FUNS(cvec, cmat, zeros_c, std::complex<double>(0.0))
CREATE_EYE_FUN(mat, eye, 0.0, 1.0)
CREATE_EYE_FUN(bmat, eye_b, bin(0), bin(1))
CREATE_EYE_FUN(imat, eye_i, 0, 1)
CREATE_EYE_FUN(cmat, eye_c, std::complex<double>(0.0), std::complex<double>(1.0))
vec impulse(int size)
{
vec t(size);
t.clear();
t[0] = 1.0;
return t;
}
vec linspace(double from, double to, int points)
{
if (points < 2) {
// This is the "Matlab definition" of linspace
vec output(1);
output(0) = to;
return output;
}
else {
vec output(points);
double step = (to - from) / double(points - 1);
int i;
for (i = 0; i < (points-1); i++)
output(i) = from + i * step;
output(i) = to;
return output;
}
}
vec zigzag_space(double t0, double t1, int K)
{
it_assert(K > 0, "zigzag_space:() K must be positive");
ivec N = "0 1";
int n = 2;
for (int k = 0; k < K; k++) {
ivec Nn = 2 * N;
for (int i = 1; i < length(Nn); i += 2) {
Nn = concat(Nn, i);
n++;
}
N = Nn;
}
vec T0 = linspace(t0, t1, n);
vec Tt = zeros(n);
for (int i = 0; i < n; i++) {
Tt(i) = T0(N(i));
}
return Tt;
}
// Construct a Hadamard-imat of size "size"
imat hadamard(int size)
{
it_assert(size > 0, "hadamard(): size is not a power of 2");
int logsize = ceil_i(::log2(static_cast<double>(size)));
it_assert(pow2i(logsize) == size, "hadamard(): size is not a power of 2");
imat H(size, size);
H(0, 0) = 1;
for (int i = 0; i < logsize; ++i) {
int pow2 = 1 << i;
for (int k = 0; k < pow2; ++k) {
for (int l = 0; l < pow2; ++l) {
H(k, l) = H(k, l);
H(k + pow2, l) = H(k, l);
H(k, l + pow2) = H(k, l);
H(k + pow2, l + pow2) = (-1) * H(k, l);
}
}
}
return H;
}
imat jacobsthal(int p)
{
int quadratic_residue;
imat out(p, p);
int i, j;
out = -1; // start with all elements equal to "-1"
// Generate a complete list of quadratic residues
for (i = 0; i < (p - 1) / 2; i++) {
quadratic_residue = ((i + 1) * (i + 1)) % p;
// set this element in all rows (col-row) = quadratic_residue
for (j = 0; j < p; j++) {
out(j, (j + quadratic_residue) % p) = 1;
}
}
// set diagonal elements to zero
for (i = 0; i < p; i++) {
out(i, i) = 0;
}
return out;
}
imat conference(int n)
{
it_assert_debug(n % 4 == 2, "conference(int n); wrong size");
int pm = n - 1; // p must be odd prime, not checked
imat out(n, n);
out.set_submatrix(1, 1, jacobsthal(pm));
out.set_submatrix(0, 0, 1, n - 1, 1);
out.set_submatrix(1, n - 1, 0, 0, 1);
out(0, 0) = 0;
return out;
}
const cmat toeplitz(const cvec &c)
{
int s = c.size();
cmat output(s, s);
cvec c_conj = conj(c);
for (int i = 1; i < s; ++i) {
for (int j = 0; j < s - i; ++j) {
output(i + j, j) = c_conj(i);
}
}
// start from j = 0 here, because the main diagonal is not conjugated
for (int j = 0; j < s; ++j) {
for (int i = 0; i < s - j; ++i) {
output(i, i + j) = c(j);
}
}
return output;
}
mat rotation_matrix(int dim, int plane1, int plane2, double angle)
{
mat m;
double c = std::cos(angle), s = std::sin(angle);
it_assert(plane1 >= 0 && plane2 >= 0 &&
plane1 < dim && plane2 < dim && plane1 != plane2,
"Invalid arguments to rotation_matrix()");
m.set_size(dim, dim, false);
m = 0.0;
for (int i = 0; i < dim; i++)
m(i, i) = 1.0;
m(plane1, plane1) = c;
m(plane1, plane2) = -s;
m(plane2, plane1) = s;
m(plane2, plane2) = c;
return m;
}
void house(const vec &x, vec &v, double &beta)
{
double sigma, mu;
int n = x.size();
v = x;
if (n == 1) {
v(0) = 1.0;
beta = 0.0;
return;
}
sigma = sum(sqr(x(1, n - 1)));
v(0) = 1.0;
if (sigma == 0.0)
beta = 0.0;
else {
mu = std::sqrt(sqr(x(0)) + sigma);
if (x(0) <= 0.0)
v(0) = x(0) - mu;
else
v(0) = -sigma / (x(0) + mu);
beta = 2 * sqr(v(0)) / (sigma + sqr(v(0)));
v /= v(0);
}
}
void givens(double a, double b, double &c, double &s)
{
double t;
if (b == 0) {
c = 1.0;
s = 0.0;
}
else {
if (fabs(b) > fabs(a)) {
t = -a / b;
s = -1.0 / std::sqrt(1 + t * t);
c = s * t;
}
else {
t = -b / a;
c = 1.0 / std::sqrt(1 + t * t);
s = c * t;
}
}
}
void givens(double a, double b, mat &m)
{
double t, c, s;
m.set_size(2, 2);
if (b == 0) {
m(0, 0) = 1.0;
m(1, 1) = 1.0;
m(1, 0) = 0.0;
m(0, 1) = 0.0;
}
else {
if (fabs(b) > fabs(a)) {
t = -a / b;
s = -1.0 / std::sqrt(1 + t * t);
c = s * t;
}
else {
t = -b / a;
c = 1.0 / std::sqrt(1 + t * t);
s = c * t;
}
m(0, 0) = c;
m(1, 1) = c;
m(0, 1) = s;
m(1, 0) = -s;
}
}
mat givens(double a, double b)
{
mat m(2, 2);
givens(a, b, m);
return m;
}
void givens_t(double a, double b, mat &m)
{
double t, c, s;
m.set_size(2, 2);
if (b == 0) {
m(0, 0) = 1.0;
m(1, 1) = 1.0;
m(1, 0) = 0.0;
m(0, 1) = 0.0;
}
else {
if (fabs(b) > fabs(a)) {
t = -a / b;
s = -1.0 / std::sqrt(1 + t * t);
c = s * t;
}
else {
t = -b / a;
c = 1.0 / std::sqrt(1 + t * t);
s = c * t;
}
m(0, 0) = c;
m(1, 1) = c;
m(0, 1) = -s;
m(1, 0) = s;
}
}
mat givens_t(double a, double b)
{
mat m(2, 2);
givens_t(a, b, m);
return m;
}
template void eye(int, mat &);
template void eye(int, bmat &);
template void eye(int, imat &);
template void eye(int, cmat &);
template vec linspace_fixed_step(double, double, double);
template ivec linspace_fixed_step(int, int, int);
template svec linspace_fixed_step(short int, short int, short int);
} // namespace itpp
matfunc.h
Various functions on vectors and matrices - header file.
itpp::zeros_i
ITPP_EXPORT ivec zeros_i(int size)
A Int vector of zeros.
itpp::ones_b
ITPP_EXPORT bvec ones_b(int size)
A Binary vector of ones.
itpp::ones_i
ITPP_EXPORT ivec ones_i(int size)
A Int vector of ones.
itpp::eye_c
ITPP_EXPORT cmat eye_c(int size)
A Double Complex (size,size) unit matrix.
itpp::ceil_i
int ceil_i(double x)
The nearest larger integer.
Definition: converters.h:339
itpp::impulse
vec impulse(int size)
Impulse vector.
Definition: specmat.cpp:98
itpp::givens_t
void givens_t(double a, double b, mat &m)
Calculate the transposed Givens rotation matrix.
Definition: specmat.cpp:336
itpp::find
ivec find(const bvec &invector)
Return a integer vector with indicies where bvec == 1.
Definition: specmat.cpp:40
itpp::conj
cvec conj(const cvec &x)
Conjugate of complex value.
Definition: elem_math.h:226
itpp::zeros_b
ITPP_EXPORT bvec zeros_b(int size)
A Binary vector of zeros.
itpp::eye
template void eye(int, mat &)
Template instantiation of eye.
itpp::sum
T sum(const Vec< T > &v)
Sum of all elements in the vector.
Definition: matfunc.h:59
it_assert
#define it_assert(t, s)
Abort if t is not true.
Definition: itassert.h:94
log_exp.h
Logarithmic and exponenential functions - header file.
itpp::length
int length(const Vec< T > &v)
Length of vector.
Definition: matfunc.h:51
it_assert_debug
#define it_assert_debug(t, s)
Abort if t is not true and NDEBUG is not defined.
Definition: itassert.h:107
itpp::zeros_c
ITPP_EXPORT cvec zeros_c(int size)
A Double Complex vector of zeros.
itpp::toeplitz
const cmat toeplitz(const cvec &c)
Generate symmetric Toeplitz matrix from vector c (complex valued)
Definition: specmat.cpp:210
itpp::sin
vec sin(const vec &x)
Sine function.
Definition: trig_hyp.h:54
itpp::zeros
ITPP_EXPORT vec zeros(int size)
A Double vector of zeros.
itpp::angle
vec angle(const cvec &x)
Angle.
Definition: elem_math.h:218
itpp::log2
vec log2(const vec &x)
log-2 of the elements
Definition: log_exp.cpp:36
itpp::pow2i
int pow2i(int x)
Calculate two to the power of x (2^x); x is integer.
Definition: log_exp.h:53
itpp::pow2
double pow2(double x)
Calculate two to the power of x (2^x)
Definition: log_exp.h:55
itpp::size
int size(const Vec< T > &v)
Length of vector.
Definition: matfunc.h:55
itpp
itpp namespace
Definition: itmex.h:36
itpp::conference
imat conference(int n)
Conference matrix.
Definition: specmat.cpp:196
itpp::jacobsthal
imat jacobsthal(int p)
Jacobsthal matrix.
Definition: specmat.cpp:172
itpp::ones
ITPP_EXPORT vec ones(int size)
A float vector of ones.
itpp::rotation_matrix
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
itpp::house
void house(const vec &x, vec &v, double &beta)
Calcualte the Householder vector.
Definition: specmat.cpp:251
itcompat.h
IT++ compatibility types and functions.
itpp::sqr
vec sqr(const cvec &data)
Absolute square of elements.
Definition: elem_math.cpp:36
specmat.h
Definitions of special vectors and matrices.
itpp::linspace_fixed_step
template vec linspace_fixed_step(double, double, double)
Template instantiation of linspace_fixed_step.
itpp::linspace
vec linspace(double from, double to, int points)
linspace (works in the same way as the MATLAB version)
Definition: specmat.cpp:106
itpp::bin
Binary arithmetic (boolean) class.
Definition: binary.h:56
itpp::sqrt
vec sqrt(const vec &x)
Square root of the elements.
Definition: elem_math.h:123
itpp::to
T to(double x)
Convert double to T.
Definition: fix_functions.h:365
itpp::zigzag_space
vec zigzag_space(double t0, double t1, int K)
Zig-zag space function (variation on linspace)
Definition: specmat.cpp:125
itpp::eye_b
ITPP_EXPORT bmat eye_b(int size)
A Binary (size,size) unit matrix.
itpp::cos
vec cos(const vec &x)
Cosine function.
Definition: trig_hyp.h:58
itpp::ones_c
ITPP_EXPORT cvec ones_c(int size)
A float Complex vector of ones.
itpp::eye_i
ITPP_EXPORT imat eye_i(int size)
A Int (size,size) unit matrix.
elem_math.h
Elementary mathematical functions - header file.
itpp::givens
void givens(double a, double b, double &c, double &s)
Calculate the Givens rotation values.
Definition: specmat.cpp:277
itpp::hadamard
imat hadamard(int size)
Hadamard matrix.
Definition: specmat.cpp:149
bmat
Mat< bin > bmat
bin matrix
Definition: mat.h:508
itpp::concat
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

Generated on Sat Sep 30 2017 07:04:09 for IT++ by Doxygen 1.8.13