ProteoWizard
mexhat.hpp
Go to the documentation of this file.
1 //
2 // $Id: mexhat.hpp 5313 2013-12-17 18:06:54Z chambm $
3 //
4 //
5 // Original author: Witold Wolski <wewolski@gmail.com>
6 //
7 // Copyright : ETH Zurich
8 //
9 // Licensed under the Apache License, Version 2.0 (the "License");
10 // you may not use this file except in compliance with the License.
11 // You may obtain a copy of the License at
12 //
13 // http://www.apache.org/licenses/LICENSE-2.0
14 //
15 // Unless required by applicable law or agreed to in writing, software
16 // distributed under the License is distributed on an "AS IS" BASIS,
17 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18 // See the License for the specific language governing permissions and
19 // limitations under the License.
20 //
21 
22 #ifndef FILTERTYPES_H
23 #define FILTERTYPES_H
24 
25 namespace ralab{
26  namespace base{
27  namespace filter{
28  namespace utilities{
29  /*! \brief Mexican hat wavelet.
30 
31  \f{eqnarray*}{
32  T_1 &= {1 \over {\sqrt {2\pi}\sigma^3}}\\
33  T_2 &= \left( 1 - {t^2 \over \sigma^2} \right)\\
34  T_3 &= e^{-t^2 \over 2\sigma^2} \\
35  \psi(t) &= T_1 \cdot T_2 \cdot T_3
36  \f}
37  is the negative normalized second derivative of a Gaussian function. (Source Wikipedia.)
38  (LoG) Laplace of Gaussian, with Laplace's operator having a (-1,2,-1) Kernel.
39 
40  \ingroup FILTER
41  */
42  template<typename TReal >
43  struct Mexican_Hat: std::unary_function< TReal, TReal> {
44 
46  TReal mu, //!< mean
47  TReal sigma //!<standard deviation
48  )
49  :mu_(mu),
50  sigma_(sigma)
51  {}
52  /*! \brief operator */
53  TReal operator()( TReal x )
54  {
55  TReal two = TReal(2);
56  TReal t1( 1 / (pow(TReal(2. * ralab::constants::PI) , TReal(.5)) * pow(sigma_,TReal(3.))) );
57  TReal t2(1 - pow( x-mu_ , two)/ pow(sigma_, two ) );
58  TReal t3( exp(-pow((x-mu_) , two )/( 2 * pow( sigma_, two ) ) ) );
59  return( t1 * t2 * t3 );
60  }
61  protected:
62  TReal mu_;
63  TReal sigma_;
64  };
65 
66  /*! \brief Mexican hat wavelet Version 2.
67 
68  For Gaussian of the same amplitude and matching withs, generates response of same amplitude.
69 
70  \f{eqnarray*}{
71  T_1 &= {1 \over {\sigma}}\\
72  T_2 &= \left( 1 - {t^2 \over \sigma^2} \right)\\
73  T_3 &= e^{-t^2 \over 2\sigma^2}\\
74  \psi(t) &= T_1 \cdot T_2 \cdot T_3
75  \f}
76 
77  Note, the change in the Term \f$T_1\f$ compared with the Mexican_Hat functor.
78  */
79 
80  template<typename TReal >
81  struct Mexican_Hat2 : std::unary_function<TReal, TReal >{
82 
84  TReal mu, //!< mean
85  TReal sigma //!<standard deviation
86  )
87  :mu_(mu),
88  sigma_(sigma)
89  {}
90 
91  TReal operator()( TReal x )
92  {
93  TReal two = TReal(2);
94  TReal t1( 1 / sigma_ );
95  TReal t2(1 - pow( x-mu_ , two)/ pow(sigma_,two ) );
96  TReal t3( exp(-pow((x-mu_) , two )/( 2 * pow( sigma_, two ) ) ) );
97  return( t1 * t2 * t3 );
98  }
99  protected:
100  TReal mu_;
101  TReal sigma_;
102  };
103 
104 
105 
106  /*!\brief Scales a mother wavelet so that the conditions hold:
107  \f[ sum(mh) = 0 \f]
108  \f[ sum(mh^2) = 1 \f]
109  */
110  template<typename TReal>
112  std::vector<TReal> &mh, //the wavelet to scale
113  std::vector<TReal> &x //temporary vector....
114  )
115  {
116  //do this so that the sum of wavelet equals zero
117  TReal sum = std::accumulate(mh.begin() , mh.end() , 0.);
118  sum /= mh.size();
119  std::transform(mh.begin(),mh.end(),mh.begin(),std::bind2nd(std::minus<TReal>(), sum )) ;
120 
121  //compute sum of square...
122  std::transform( mh.begin(), mh.end(), x.begin(), ralab::base::stats::NthPower<2,TReal>() ); //first sqaure all elements
123  TReal sumsq = sqrt(std::accumulate(x.begin(), x.end() , TReal(0.)));
124  std::transform(mh.begin() , mh.end() , mh.begin() , std::bind2nd(std::divides<TReal>(), sumsq ) ) ;
125 
126  //this is just to verify the result
127  std::transform( mh.begin(), mh.end(), x.begin(), ralab::base::stats::NthPower<2,TReal>() ); //first sqaure all elements
128  sumsq = std::accumulate(x.begin(), x.end() , TReal(0.));
129  return sumsq;
130 
131  }
132 
133  template<typename TReal>
134  TReal getMaxHatWorker( TReal sigma, std::vector<TReal> &mh, std::vector<TReal> &x )
135  {
136  ralab::base::filter::utilities::Mexican_Hat<TReal> mexHatGenerator(0.,sigma);
137  mh.resize( x.size() );
138  std::transform( x.begin() ,x.end(),mh.begin(),mexHatGenerator);
139  scaleWavelet(mh, x);
140  TReal sum = std::accumulate(mh.begin(),mh.end(),0.);
141  return sum;
142  }
143  }//utilities
144  }//filter
145  }//base
146 }//ralab
147 
148 #endif
TReal operator()(TReal x)
operator
Definition: mexhat.hpp:53
TReal getMaxHatWorker(TReal sigma, std::vector< TReal > &mh, std::vector< TReal > &x)
Definition: mexhat.hpp:134
Mexican_Hat(TReal mu, TReal sigma)
Definition: mexhat.hpp:45
TReal scaleWavelet(std::vector< TReal > &mh, std::vector< TReal > &x)
Scales a mother wavelet so that the conditions hold: .
Definition: mexhat.hpp:111
const double PI(3.14159265358979323846264338327950288)
the ratio of the circumference of a circle to its diameter;
EQUISPACEINTERPOL Interpolation on a equidistantly spaced grid.
Definition: base.hpp:39
KernelTraitsBase< Kernel >::space_type::abscissa_type x
Mexican hat wavelet Version 2.
Definition: mexhat.hpp:81
void filter(const TContainer &data, const TContainer &filter, TContainer &result, bool circular=false, uint32_t sides=2)
Applies linear convolution (filtering) to a univariate time series.
Definition: filter.hpp:112