24 #ifndef _MATCHEDFILTER_HPP_ 25 #define _MATCHEDFILTER_HPP_ 39 namespace MatchedFilter {
42 template <
typename X,
typename Y>
54 template <
typename space_type>
67 return domain.second - domain.first;
70 abscissa_type
dx()
const 72 return samples.empty() ? 0 : domainWidth()/(samples.size()-1);
75 abscissa_type
x(
typename samples_type::size_type index)
const 77 return domain.first + domainWidth() * index / (samples.size()-1);
80 typename samples_type::size_type
sampleIndex(
const abscissa_type&
x)
const 82 typename samples_type::size_type sampleCount = samples.size();
83 int result = (int)round((sampleCount-1)*(x-domain.first)/domainWidth());
84 if (result < 0) result = 0;
85 if (result > (
int)(sampleCount-1)) result = sampleCount-1;
86 return (
typename samples_type::size_type)(result);
89 const ordinate_type&
sample(abscissa_type
x)
const 91 return samples[sampleIndex(x)];
96 template <
typename space_type>
97 std::ostream& operator<<(std::ostream& os, const SampledData<space_type>& data)
99 os <<
"[" << data.domain.first <<
"," << data.domain.second <<
"] " 100 <<
"(" << data.samples.size() <<
" samples)\n";
103 for (
unsigned int index=0; index!=data.samples.size(); ++index, ++it)
104 os << data.x(index) <<
"\t" << *it << std::endl;
110 template <
typename Kernel>
125 template <
typename X,
typename Y>
134 template <
typename Kernel>
142 typename KernelTraitsBase<Kernel>::space_type::abscissa_type
x;
143 typename KernelTraitsBase<Kernel>::space_type::ordinate_type
y;
149 template <
typename Kernel>
153 void (KernelConcept<Kernel>::*dummy)() = &KernelConcept<Kernel>::check;
158 template <
typename Y>
166 : dot(_dot), e2(_e2), tan2angle(_tan2angle)
169 double angle()
const {
return atan(sqrt(tan2angle))*180/M_PI;}
174 std::ostream& operator<<(std::ostream& os, const Correlation<Y>& c)
176 os <<
"<" << c.dot <<
", " << c.e2 <<
", " << c.angle() <<
">";
181 template <
typename Kernel>
204 template <
typename Kernel>
211 checkKernelConcept<Kernel>();
214 for (
int i=-sampleRadius; i<=sampleRadius; i++)
215 filter.push_back(kernel(i*dx - shift));
221 inline double norm(
double d) {
return d*d;}
222 inline double conj(
double d) {
return d;}
225 template <
typename Filter>
228 double normalization = 0;
229 for (
typename Filter::const_iterator it=filter.begin(); it!=filter.end(); ++it)
230 normalization +=
norm(*it);
231 normalization = sqrt(normalization);
233 for (
typename Filter::iterator it=filter.begin(); it!=filter.end(); ++it)
234 *it /= normalization;
238 template <
typename Kernel>
239 std::vector<typename KernelTraits<Kernel>::filter_type>
245 checkKernelConcept<Kernel>();
248 std::vector<filter_type> filters;
250 for (
int i=0; i<subsampleFactor; i++)
251 filters.push_back(
createFilter(kernel, sampleRadius, dx, dx*i/subsampleFactor));
253 for_each(filters.begin(), filters.end(), normalizeFilter<filter_type>);
259 template<
typename Kernel>
266 checkKernelConcept<Kernel>();
270 for (; samples!=samplesEnd; ++samples, ++
filter)
272 result.
dot += (*samples) *
conj(*filter);
273 normData +=
norm(*samples);
276 double normDot =
norm(result.
dot);
277 result.
e2 = (std::max)(normData - normDot, 0.);
278 result.
tan2angle = normDot>0 ? result.
e2/normDot : std::numeric_limits<double>::infinity();
285 template<
typename Kernel>
288 const Kernel& kernel,
292 checkKernelConcept<Kernel>();
298 if (data.
samples.empty())
return result;
299 result.samples.resize((data.
samples.size()-1) * subsampleFactor + 1);
309 unsigned int sampleIndex = sampleRadius;
310 for (
typename samples_type::const_iterator itData = data.
samples.begin() + sampleRadius;
311 itData + sampleRadius != data.
samples.end(); ++itData, ++sampleIndex)
312 for (
unsigned int filterIndex=0; filterIndex<filters.size(); ++filterIndex)
314 unsigned int index = sampleIndex * filters.size() + filterIndex;
316 if (index >= result.samples.size())
319 details::computeCorrelation<Kernel>(itData-sampleRadius, itData+sampleRadius+1,
320 filters[filterIndex].begin(),
321 result.samples[index]);
333 #endif // _MATCHEDFILTER_HPP_
ProductSpace< double, double > DxD
Kernel::space_type space_type
void normalizeFilter(Filter &filter)
ProductSpace< X, Y > space_type
const ordinate_type & sample(abscissa_type x) const
std::pair< abscissa_type, abscissa_type > domain_type
ProductSpace< double, std::complex< double > > DxCD
KernelTraits< Kernel >::filter_type createFilter(const Kernel &kernel, int sampleRadius, typename KernelTraits< Kernel >::abscissa_type dx, typename KernelTraits< Kernel >::abscissa_type shift)
SampledData< correlation_space_type > correlation_data_type
void computeCorrelation(typename KernelTraits< Kernel >::samples_type::const_iterator samples, typename KernelTraits< Kernel >::samples_type::const_iterator samplesEnd, typename KernelTraits< Kernel >::samples_type::const_iterator filter, typename KernelTraits< Kernel >::correlation_type &result)
std::vector< ordinate_type > samples_type
SampledData< space_type > sampled_data_type
samples_type::size_type sampleIndex(const abscissa_type &x) const
abscissa_type domainWidth() const
void checkKernelConcept()
KernelTraits< Kernel >::correlation_data_type computeCorrelationData(const typename KernelTraits< Kernel >::sampled_data_type &data, const Kernel &kernel, int sampleRadius, int subsampleFactor)
std::vector< typename KernelTraits< Kernel >::filter_type > createFilters(const Kernel &kernel, int sampleRadius, int subsampleFactor, typename KernelTraits< Kernel >::abscissa_type dx)
sampled_data_type::samples_type samples_type
KernelTraitsBase< Kernel >::space_type space_type
space_type::ordinate_type ordinate_type
Correlation(Y _dot=0, double _e2=0, double _tan2angle=0)
space_type::abscissa_type abscissa_type
Correlation< ordinate_type > correlation_type
space_type::ordinate_type ordinate_type
space_type::abscissa_type abscissa_type
Dummy< &checkKernelConcept< Kernel > > dummy
KernelTraitsBase< Kernel >::space_type::abscissa_type x
KernelTraitsBase< Kernel >::space_type::ordinate_type y
ProductSpace< abscissa_type, correlation_type > correlation_space_type
abscissa_type x(typename samples_type::size_type index) const
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.