SigmaTransform
The SigmaTransform unifies various known signal processing transforms, like the STFT and the wavelet transform, differing only by a specific diffeomorphism.
 All Classes Namespaces Files Functions Variables Typedefs Friends Macros
SigmaTransform_util.h
Go to the documentation of this file.
1 #ifndef SIGMATRANSFORM_UTIL_H
2 #define SIGMATRANSFORM_UTIL_H
3 
4 #include <vector>
5 #include <map>
6 #include <complex>
7 #include <thread>
8 #include <array>
9 #include <iostream>
10 #include <fstream>
11 #include <algorithm>
12 #include <mutex>
13 #include <functional>
14 
15 #define DEBUG
16 
17 #ifdef DEBUG
18 
19 #include <chrono>
20 #include <iomanip>
21 
22 #endif
23 
24 #include <fftw3.h>
25 
26 namespace SigmaTransform {
27 
28  const double PI = 3.141592654;
29  const double pow2_1_4 = 1.189207115;
30 
31  #ifdef DEBUG
32 
35  class Chronometer {
36  std::chrono::steady_clock::time_point _tic;
37  public:
38  Chronometer();
39 
44  void tic();
45 
52  Chronometer& toc(std::string const& str );
53  };
54 
55  #endif
56 
66  template<size_t N>
67  class point {
68  // private, internally handled data
69  std::array<double,N> _data;
70 
71  public:
72 
73  // constructors
74  point<N>() : point<N>({0}) {}
75  point<N>(double const& x) { for( auto& _x : _data ) _x = x; }
76  point<N>( std::array<double,N> dat ) : _data( dat ) {}
77 
78  // we may overwrite a point
79  void operator=( std::array<int,N> dat ) { for(int k=0;k<N;++k)_data[k]=(double)dat[k]; }
80 
81  // overloaded arithmetic operator with const-qualifier
82  point<N> operator+( point<N> const& r ) const { point<N> p(_data); auto R=r.begin(); for( auto& x : p ) x += *R++; return std::move( p ); }
83  point<N> operator-( point<N> const& r ) const { point<N> p(_data); auto R=r.begin(); for( auto& x : p ) x -= *R++; return std::move( p ); }
84  point<N> operator*( point<N> const& r ) const { point<N> p(_data); auto R=r.begin(); for( auto& x : p ) x *= *R++; return std::move( p ); }
85  point<N> operator/( point<N> const& r ) const { point<N> p(_data); auto R=r.begin(); for( auto& x : p ) x /= *R++; return std::move( p ); }
86  point<N> operator/( double const& r ) const { point<N> p(_data); for( auto& x : p ) x /= r; return std::move( p ); }
87  point<N> operator*( double const& r ) const { point<N> p(_data); for( auto& x : p ) x *= r; return std::move( p ); }
88  point<N> operator+( double const& r ) const { point<N> p(_data); for( auto& x : p ) x += r; return std::move( p ); }
89  point<N> operator-( double const& r ) const { point<N> p(_data); for( auto& x : p ) x -= r; return std::move( p ); }
90 
91  // overloaded arithmetic operator
92  void operator+=( point<N> const& r ) { auto R=r.begin(); for( auto& x : _data ) x += *R++; }
93  void operator-=( point<N> const& r ) { auto R=r.begin(); for( auto& x : _data ) x -= *R++; }
94  void operator*=( point<N> const& r ) { auto R=r.begin(); for( auto& x : _data ) x *= *R++; }
95  void operator/=( point<N> const& r ) { auto R=r.begin(); for( auto& x : _data ) x /= *R++; }
96 
97  // sum / prod / abs/ sq(uare)
98  double sum() const { double s = 0; for( auto const& x : _data ) s+=x; return s; }
99  double prod() const { double p = 1; for( auto const& x : _data ) p*=x; return p; }
100  point<N> sq() const { point<N> p( _data ); for( auto& x : p ) x *= x; return std::move( p ); }
101  point<N> abs() const { point<N> p( _data ); for( auto& x : p ) x = (x>=0)?x:-x; return std::move( p ); }
102 
103  // apply specific functions to each component in the "point-vector"
104  point<N> apply( double(*fun)(double const&) ) const { point<N> p; for(int k=0;k<N;++k) p[k] = fun(_data[k]); return std::move( p ); }
105  point<N> apply( double(*fun)(double) ) const { point<N> p; for(int k=0;k<N;++k) p[k] = fun(_data[k]); return std::move( p ); }
106 
107  // comparison-related operator
108  point<N> operator>( double const& d ) const { point<N> p(_data); for( auto& x : p ) x = x>d; return std::move( p ); }
109  point<N> operator<( double const& d ) const { point<N> p(_data); for( auto& x : p ) x = x<d; return std::move( p ); }
110  point<N> operator<=( double const& d ) const { point<N> p(_data); for( auto& x : p ) x = x<=d; return std::move( p ); }
111  point<N> operator>=( double const& d ) const { point<N> p(_data); for( auto& x : p ) x = x>=d; return std::move( p ); }
112  point<N> operator==( double const& d ) const { point<N> p(_data); for( auto& x : p ) x = x==d; return std::move( p ); }
113  point<N> operator>( point<N> const& r ) const { point<N> p(_data); auto R=r.begin(); for(auto&x : p) x = x>*R++; return std::move( p ); }
114  point<N> operator<( point<N> const& r ) const { point<N> p(_data); auto R=r.begin(); for(auto&x : p) x = x<*R++; return std::move( p ); }
115  point<N> operator>=( point<N> const& r ) const { point<N> p(_data); auto R=r.begin(); for(auto&x : p) x = x>=*R++; return std::move( p ); }
116  point<N> operator<=( point<N> const& r ) const { point<N> p(_data); auto R=r.begin(); for(auto&x : p) x = x<=*R++; return std::move( p ); }
117  point<N> operator==( point<N> const& r ) const { point<N> p(_data); auto R=r.begin(); for(auto&x : p) x = x==*R++; return std::move( p ); }
118  operator bool() { bool b=true; for(auto&x:_data) if(x!=1) b=false; return b; }
119 
120  // acquire internal data
121  double& operator[](size_t const& ind) { return _data[ind]; }
122  double const& operator[](size_t const& ind) const { return _data[ind]; }
123  double* get() { return _data; }
124 
125  // iterators, for compatibility with the C++ - Std Template Library
126  typedef typename std::array<double,N>::iterator iterator;
127  typedef typename std::array<double,N>::const_iterator const_iterator;
128  typedef typename std::array<double,N>::reverse_iterator reverse_iterator;
129  typedef typename std::array<double,N>::const_reverse_iterator const_reverse_iterator;
130 
131  iterator begin() { return _data.begin(); }
132  iterator end() { return _data.end(); }
133 
134  const_iterator begin() const { return _data.begin(); }
135  const_iterator end() const { return _data.end(); }
136 
137  const_iterator cbegin() const { return _data.cbegin(); }
138  const_iterator cend() const { return _data.cend(); }
139 
140  reverse_iterator rbegin() { return _data.rbegin(); }
141  reverse_iterator rend() { return _data.rend(); }
142 
143  const_reverse_iterator crbegin() const { return _data.crbegin(); }
144  const_reverse_iterator crend() const { return _data.crend(); }
145 
146  // friends for streams
147  friend std::ostream& operator<<( std::ostream& os , point<N> const& p ) { os<<"("; for( auto& x : p ) os << x << " "; return os << "\b)"; }
148  friend std::istream& operator>>( std::istream& is , point<N>& p ) { for( auto& x : p ) is >> x; return is; }
149 
150  };
151 
152  // shorthands for cleaner code
153  using cmpx = std::complex<double>;
154  using cxVec = std::vector<cmpx>;
155 
156  template<size_t N> using diffFunc = std::function<point<N>(const point<N>&)>;
157  template<size_t N> using winFunc = std::function<cmpx(const point<N>&)>;
158  template<size_t N> using actFunc = std::function<point<N>(const point<N>&,const point<N>&)>;
159 
168  cxVec loadAscii1D( std::string const& filename );
169 
180  cxVec loadAscii2D( std::string const& filename, int &x , int &y );
181 
190  std::vector<double> linspace( const double &L , const double &R , const int &N );
191 
199  std::vector<double> FourierAxis( const double &fs , const unsigned &len );
200 
208  void StartRecursiveLoop( const std::vector<int>& max , std::function<void(const std::vector<int>&)> toDo );
209 
220  void recursiveLoop( std::vector<int> &index, const std::vector<int>& max , const int& curr_ind ,
221  std::function<void(const std::vector<int>&)> toDo );
222 
229  template<size_t N>
230  std::vector<point<N>> meshgridN( std::array<std::vector<double>,N> const& dom ) {
231  // make sizeVector for loop
232  std::vector<int> sizeVec(0);
233  int p = 1;
234  for( auto& d : dom ) {
235  p *= d.size();
236  sizeVec.push_back( d.size() );
237  }
238  // reserve space
239  std::vector<point<N>> out; out.reserve( p );
240  // use recursive loop - very slow, but replaces up to n nested loops
241  StartRecursiveLoop( sizeVec , [&](std::vector<int> const& index ) {
242  // get current point
243  point<N> p;
244  for( int k = 0 ; k < index.size() ; ++k )
245  p[k] = dom[k][index[k]];
246  // insert into vectors
247  out.emplace_back( p );
248  } );
249  // return output vector
250  return std::move( out );
251  }
252 
259  template<size_t N>
260  std::vector<point<N>> meshgridN( std::vector<double> const& dom ) {
261  // make sizeVector for loop
262  std::vector<int> sizeVec(0);
263  int p = 1;
264  for(int k = 0 ; k < N;++k) {
265  p *= dom.size();
266  sizeVec.push_back( dom.size() );
267  }
268  // reserve space
269  std::vector<point<N>> out; out.reserve( p );
270  // use recursive loop - very slow, but replaces up to n nested loops
271  StartRecursiveLoop( sizeVec , [&](std::vector<int> const& index ) {
272  // get current point
273  point<N> p;
274  for( int k = 0 ; k < index.size() ; ++k )
275  p[k] = dom[index[k]];
276  // insert into vectors
277  out.emplace_back( p );
278  } );
279  // return output vector
280  return std::move( out );
281  }
282 
291  template<size_t N>
292  void save2file_asc( std::string const& filename , const cxVec &out , const point<N> & dim = point<N>(0) ) {
293  // use a stringstream buffer
294  std::stringstream ss;
295  // set precision
296  ss.precision(6);
297  //write to buffer
298  for( const auto &c : out )
299  ss << c.real() << "," << c.imag() << std::endl;
300  // output filestream
301  std::ofstream os(filename);
302  // save dimensions of vector
303  if( dim.prod() == out.size() )
304  os << "dim=" << dim << "\n";
305  else
306  os << "dim=(" <<out.size()<<")\n";
307  // write buffer to file
308  os<<ss.str();
309  }
310 
318  void save2file_bin( std::string const& filename , const cxVec &out );
319 
327  template<size_t N>
328  point<N> plus( const point<N> &l , const point<N> &r ) { return l + r; }
329 
337  template<size_t N>
338  point<N> minus( const point<N> &l , const point<N> &r ) { return l - r; }
339 
346  template<size_t N>
347  point<N> logabs( const point<N> &x ) { return (x.abs()+1E-16).apply( log2 ); }
348 
355  template<size_t N>
356  point<N> logpos( const point<N> &x ) { return (x>0).apply( log2 ); }
357 
364  template<size_t N>
365  point<N> id( const point<N> &x ) { return x; }
366 
373  template<size_t N>
374  cmpx gauss( const point<N> &x ) {
375  return exp( -PI * x.sq().sum() ) * point<N>(pow2_1_4).prod();
376  }
377 
385  template<size_t N>
386  cmpx gauss_stddev( const point<N> &x , const point<N> &stddev ) {
387  return exp( -PI * (x/stddev).sq().sum() ) * point<N>(pow2_1_4).prod();
388  }
389 
390 } // namespace SigmaTransform
391 
392 #endif //SIGMATRANSFORM_H