ProteoWizard
mexhat.hpp
Go to the documentation of this file.
1//
2// $Id: $
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
25namespace 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 {
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
KernelTraitsBase< Kernel >::space_type::abscissa_type x
TReal getMaxHatWorker(TReal sigma, std::vector< TReal > &mh, std::vector< TReal > &x)
Definition mexhat.hpp:134
TReal scaleWavelet(std::vector< TReal > &mh, std::vector< TReal > &x)
Scales a mother wavelet so that the conditions hold:
Definition mexhat.hpp:111
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
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:40
Mexican hat wavelet Version 2.
Definition mexhat.hpp:81