Vespucci  1.0.0
accessory_impl.h
Go to the documentation of this file.
1 /*******************************************************************************
2  Copyright (C) 2014-2016 Wright State University - All Rights Reserved
3  Daniel P. Foose - Maintainer/Lead Developer
4 
5  This file is part of Vespucci.
6 
7  Vespucci is free software: you can redistribute it and/or modify
8  it under the terms of the GNU General Public License as published by
9  the Free Software Foundation, either version 3 of the License, or
10  (at your option) any later version.
11 
12  Vespucci is distributed in the hope that it will be useful,
13  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15  GNU General Public License for more details.
16 
17  You should have received a copy of the GNU General Public License
18  along with Vespucci. If not, see <http://www.gnu.org/licenses/>.
19 *******************************************************************************/
20 #ifndef ACCESSORY_IMPL_H
21 #define ACCESSORY_IMPL_H
22 //Implementation of template functions
24 #include <Global/libvespucci.h>
25 namespace Vespucci
26 {
27  namespace Math
28  {
29  template <typename T> T conv_fft(const T &A, const T &B, std::string type);
30  template <typename T> T rotate(const T &X, arma::uword shift_by, bool slide_back = true);
31  template <typename T> T diff(const T &X, arma::uword deriv_order = 1);
32  }
33 }
34 
35 
47 template <typename T> T Vespucci::Math::diff(const T &X, arma::uword deriv_order)
48 {
49  T difference = X;
50  if (deriv_order > 0){ //perform first derivative of input
51  for (arma::uword i = 0; i < deriv_order; ++i){
52  T offset(difference.n_rows, difference.n_cols);
53  offset.zeros();
54  offset.rows(1, offset.n_rows - 1) = difference.rows(0, difference.n_rows - 2);
55  difference = difference - offset;
56  difference.shed_row(0); //difference will have one less row
57  }
58  }
59 
60  return difference;
61 }
62 
63 
64 
65 template <typename T> T Vespucci::Math::rotate(const T &X, arma::uword shift_by, bool slide_back)
66 {
67  T start;
68  T end;
69 
70  if (shift_by >= X.n_rows){
71  shift_by = shift_by - X.n_rows;
72  }
73 
74  if (shift_by == 0){
75  return X;
76  }
77  else if (shift_by == 1){
78  if (slide_back){
79  start = X.rows(1, X.n_rows - 1);
80  end = X.row(0);
81  }
82  else{
83  start = X.row(X.n_rows - 1);
84  end = X.rows(0, X.n_rows - 2);
85  }
86  }
87  else{
88  if (slide_back){
89  start = X.rows(shift_by, X.n_rows - 1);
90  end = X.rows(0, shift_by - 1);
91  }
92  else{
93  start = X.rows(X.n_rows - shift_by, X.n_rows - 1);
94  end = X.rows(0, X.n_rows - shift_by - 1);
95  }
96  }
97 
98  return join_vert(start, end);
99 }
100 
101 //Column-wise 1-D convolution (may re-write with 2d fft later)
102 template <typename T> T Vespucci::Math::conv_fft(const T &A, const T &B, std::string type)
103 {
104 
105  if (A.n_rows < B.n_rows){
106  throw std::invalid_argument("B must have the same number of rows or fewer rows than A");
107  }
108 
109  if (type == "filter"){
110  arma::uword next_p2 = std::ceil(std::log(A.n_elem) / std::log(2));
111  arma::uword n = std::pow(2, next_p2);
112  arma::cx_mat A_hat = fft(A, n);
113  arma::cx_mat B_hat = fft(B, n);
114  arma::cx_mat w = ifft(A_hat % B_hat, n);
115  w = w.rows(0, A.n_rows - 1);
116  return arma::real(w);
117  }
118  else if (type == "circular"){
119 
120  T A_padded = join_vert(A, arma::zeros(B.n_rows- 1));
121  T B_padded = join_vert(B, arma::zeros(A.n_rows - 1));
122 
123  arma::uword next_p2 = std::ceil(std::log(A.n_elem) / std::log(2));
124  arma::uword n = std::pow(2, next_p2);
125 
126  arma::cx_mat A_hat = fft(A_padded, n);
127  arma::cx_mat B_hat = fft(B_padded, n);
128  arma::cx_mat w = ifft(A_hat % B_hat, n);
129  return arma::real(w);
130  }
131  else
132  throw(std::invalid_argument("argument 'type' invalid"));
133 }
134 
135 
136 #endif // ACCESSORY_IMPL_H
T diff(const T &X, arma::uword deriv_order=1)
Vespucci::Math::diff.
T conv_fft(const T &A, const T &B, std::string type)
T rotate(const T &X, arma::uword shift_by, bool slide_back=true)
cmplx FADDEEVA() w(cmplx z, double relerr)
Definition: Faddeeva.cpp:680
A namespace for "global" functions, including math functions.