[med-svn] [SCM] aghermann branch, master, updated. 06bda7dfaa687aaf0708a024d192024e2cd58421
Andrei Zavada
johnhommer at gmail.com
Thu Jan 24 00:43:30 UTC 2013
The following commit has been merged in the master branch:
commit 9d5ae063248e9aa94c8fab83ce2d627f051dc04b
Author: Andrei Zavada <johnhommer at gmail.com>
Date: Wed Jan 9 19:32:18 2013 +0200
patterns refactoring (part 1/2)
diff --git a/src/common/alg.hh b/src/common/alg.hh
index 692ce22..b6f4301 100644
--- a/src/common/alg.hh
+++ b/src/common/alg.hh
@@ -113,7 +113,6 @@ struct SSpan {
-
template <typename T>
int
sign( T x)
diff --git a/src/sigproc/Makefile.am b/src/sigproc/Makefile.am
index eab6284..ec567da 100644
--- a/src/sigproc/Makefile.am
+++ b/src/sigproc/Makefile.am
@@ -8,6 +8,7 @@ liba_a_SOURCES := \
exstrom.cc exstrom.hh \
ext-filters.cc ext-filters.hh ext-filters.ii \
sigproc.cc sigproc.hh sigproc.ii \
+ patterns.cc patterns.hh patterns.ii \
winfun.cc winfun.hh
if DO_PCH
@@ -15,7 +16,8 @@ BUILT_SOURCES := \
ext-filters.hh.gch \
exstrom.hh.gch \
winfun.hh.gch \
- sigproc.hh.gch
+ sigproc.hh.gch \
+ patterns.hh.gch
%.hh.gch: %.hh
$(CXXCOMPILE) -c $<
CLEANFILES := \
diff --git a/src/sigproc/patterns.cc b/src/sigproc/patterns.cc
new file mode 100644
index 0000000..9e7f989
--- /dev/null
+++ b/src/sigproc/patterns.cc
@@ -0,0 +1,31 @@
+// ;-*-C++-*-
+/*
+ * File name: sigproc/patterns.cc
+ * Project: Aghermann
+ * Author: Andrei Zavada <johnhommer at gmail.com>
+ * Initial version: 2013-01-09
+ *
+ * Purpose: CPattern explicit pattern instantiations be here
+ *
+ * License: GPL
+ */
+
+#include "patterns.hh"
+
+using namespace std;
+
+template sigproc::CPattern<TFloat>::CPattern( const SSignalRef<TFloat>&, size_t, size_t,
+ const SPatternPPack<TFloat>&);
+template size_t sigproc::CPattern<TFloat>::find( const valarray<TFloat>&,
+ const valarray<TFloat>&,
+ const valarray<TFloat>&,
+ const valarray<TFloat>&,
+ ssize_t, int);
+template size_t sigproc::CPattern<TFloat>::find( const SSignalRef<TFloat>&,
+ ssize_t, int);
+template size_t sigproc::CPattern<TFloat>::find( const valarray<TFloat>&,
+ ssize_t, int);
+
+
+
+// eof
diff --git a/src/sigproc/patterns.hh b/src/sigproc/patterns.hh
new file mode 100644
index 0000000..b547cce
--- /dev/null
+++ b/src/sigproc/patterns.hh
@@ -0,0 +1,157 @@
+// ;-*-C++-*-
+/*
+ * File name: sigproc/patterns.hh
+ * Project: Aghermann
+ * Author: Andrei Zavada <johnhommer at gmail.com>
+ * Initial version: 2013-01-09
+ *
+ * Purpose: class CPattern
+ *
+ * License: GPL
+ */
+
+#ifndef _SIGPROC_PATTERNS_H
+#define _SIGPROC_PATTERNS_H
+
+#include <valarray>
+#include <array>
+
+#include "sigproc.hh"
+
+#if HAVE_CONFIG_H && !defined(VERSION)
+# include "config.h"
+#endif
+
+using namespace std;
+
+namespace sigproc {
+
+template <typename T>
+struct TMatch : public array<T, 4> {
+ TMatch (T _1, T _2, T _3, T _4)
+ : array<T, 4> {{_1, _2, _3, _4}}
+ {}
+ TMatch<T> () = default;
+ TMatch<T>& operator/( T dvsr)
+ {
+ return ((array<T, 4>)(*this)) / dvsr;
+ }
+ bool good_enough( const TMatch<T>& rv) const
+ {
+ for ( size_t i = 0; i < 4; ++i )
+ if ( (*this)[i] > rv[i] )
+ return false;
+ return true;
+ }
+};
+
+template <typename T>
+struct SPatternPPack {
+ int env_tightness;
+ double bwf_ffrom,
+ bwf_fupto;
+ int bwf_order;
+ double dzcdf_step,
+ dzcdf_sigma;
+ int dzcdf_smooth;
+ bool operator==( const SPatternPPack<T>& rv) const // cannot be defaulted!
+ {
+ return env_tightness == rv.env_tightness &&
+ bwf_ffrom == rv.bwf_ffrom &&
+ bwf_fupto == rv.bwf_fupto &&
+ bwf_order == rv.bwf_order &&
+ dzcdf_step == rv.dzcdf_step &&
+ dzcdf_sigma == rv.dzcdf_sigma &&
+ dzcdf_smooth == rv.dzcdf_smooth &&
+ criteria == rv.criteria;
+ }
+ TMatch<T>
+ criteria;
+}; // keep fields in order, or edit ctor by initializer_list
+
+
+
+template <typename T>
+class CPattern
+ : public SPatternPPack<T> {
+ CPattern () = delete;
+
+ public:
+ // the complete pattern signature is made of:
+ // (a) signal breadth at given tightness;
+ // (b) its course;
+ // (c) target frequency (band-passed);
+ // (d) instantaneous frequency at fine intervals;
+
+ TMatch<T>
+ match; // resulting
+
+ CPattern (const SSignalRef<T>& thing,
+ size_t ctx_before_, size_t ctx_after_,
+ const SPatternPPack<T>& Pp_)
+ : SPatternPPack<T> (Pp_),
+ match (NAN, NAN, NAN, NAN),
+ penv (thing),
+ ptarget_freq (thing),
+ pdzcdf (thing),
+ samplerate (thing.samplerate),
+ ctx_before (ctx_before_), ctx_after (ctx_after_)
+ {
+ if ( ctx_before + ctx_after >= thing.signal.size() )
+ throw invalid_argument ("pattern.size too small");
+ crit_linear_unity =
+ penv(Pp_.env_tightness).first.max() -
+ penv(Pp_.env_tightness).second.min();
+ crit_dzcdf_unity =
+ pdzcdf(Pp_.dzcdf_step,
+ Pp_.dzcdf_sigma,
+ Pp_.dzcdf_smooth).max();
+ }
+
+ size_t find( const SSignalRef<T>& field,
+ ssize_t start,
+ int inc);
+ size_t find( const valarray<T>& field,
+ ssize_t start,
+ int inc);
+ size_t find( const valarray<T>& env_u, // broken-down field
+ const valarray<T>& env_l,
+ const valarray<T>& target_freq,
+ const valarray<T>& dzcdf,
+ ssize_t start,
+ int inc);
+
+ size_t size_with_context() const
+ {
+ return ptarget_freq.signal.size();
+ }
+ size_t size_essential() const
+ {
+ return size_with_context()
+ - ctx_before - ctx_after;
+ }
+
+ private:
+ SCachedEnvelope<T>
+ penv;
+ SCachedBandPassCourse<T>
+ ptarget_freq;
+ SCachedDzcdf<T>
+ pdzcdf;
+
+ size_t samplerate;
+ size_t ctx_before,
+ ctx_after;
+
+ T crit_linear_unity;
+ double crit_dzcdf_unity;
+};
+
+#include "patterns.ii"
+
+} // namespace sigproc
+
+
+#endif
+
+// eof
diff --git a/src/sigproc/patterns.ii b/src/sigproc/patterns.ii
new file mode 100644
index 0000000..b64c1a7
--- /dev/null
+++ b/src/sigproc/patterns.ii
@@ -0,0 +1,126 @@
+// ;-*-C++-*-
+/*
+ * File name: sigproc/patterns.ii
+ * Project: Aghermann
+ * Author: Andrei Zavada <johnhommer at gmail.com>
+ * Initial version: 2013-01-09
+ *
+ * Purpose: CPattern templates
+ *
+ * License: GPL
+ */
+
+extern template CPattern<TFloat>::CPattern( const SSignalRef<TFloat>&, size_t, size_t,
+ const SPatternPPack&);
+extern template size_t CPattern<TFloat>::find( const valarray<TFloat>&,
+ const valarray<TFloat>&,
+ const valarray<TFloat>&,
+ const valarray<TFloat>&,
+ ssize_t, int);
+extern template size_t CPattern<TFloat>::find( const SSignalRef<TFloat>&,
+ ssize_t, int);
+extern template size_t CPattern<TFloat>::find( const valarray<TFloat>&,
+ ssize_t, int);
+
+
+template <typename T>
+size_t
+CPattern<T>::
+find( const valarray<T>& fenv_u,
+ const valarray<T>& fenv_l,
+ const valarray<T>& ftarget_freq,
+ const valarray<T>& fdzcdf,
+ ssize_t start,
+ int inc)
+{
+ if ( inc == 0 || inc > (ssize_t)ftarget_freq.size() ) {
+ fprintf( stderr, "CPattern::find(): bad search increment: %d\n", inc);
+ return (size_t)-1;
+ }
+
+ // printf( "course.size = %zu, fcourse.size = %zu, start = %zu\n",
+ // course.size(), fcourse.size(), start);
+ ssize_t iz = (inc > 0) ? ftarget_freq.size() - size_with_context() : 0;
+ size_t essential_part = size_essential();
+ // bool looking_further = false;
+ // T ax, bx, cx;
+ for ( ssize_t i = start; (inc > 0) ? i < iz : i > iz; i += inc ) {
+ TMatch<T>
+ diff;
+ for ( size_t j = 0; j < essential_part; ++j ) {
+ diff[0] += fdim( penv.centre( SPatternPPack<T>::env_tightness, ctx_before + j),
+ (fenv_u[i+j] + fenv_l[i+j])/2);
+ diff[1] += fdim( penv.breadth( SPatternPPack<T>::env_tightness, ctx_before + j),
+ (fenv_u[i+j] - fenv_l[i+j]));
+ diff[2] += fdim( ptarget_freq( SPatternPPack<T>::bwf_ffrom,
+ SPatternPPack<T>::bwf_fupto,
+ SPatternPPack<T>::bwf_order) [ctx_before + j],
+ ftarget_freq[i+j]);
+ diff[3] += fdim( pdzcdf( SPatternPPack<T>::dzcdf_step,
+ SPatternPPack<T>::dzcdf_sigma,
+ SPatternPPack<T>::dzcdf_smooth) [ctx_before + j],
+ fdzcdf[i+j]);
+ }
+
+ diff = diff / essential_part;
+ diff[0] /= crit_linear_unity; // normalise
+ diff[1] /= crit_linear_unity;
+ diff[2] /= crit_linear_unity;
+ diff[3] /= crit_dzcdf_unity;
+
+ // if ( i % 250 == 0 ) printf( "at %zu diff_course = %g,\tdiff_breadth = %g\t diff_dzcdf = %g\n", i, diff_course, diff_breadth, diff_dzcd);
+ if ( diff.good_enough(SPatternPPack<T>::criteria) ) {
+ // if ( !looking_further ) {
+ // looking_further = true;
+ match = diff;
+ return i;
+ }
+ }
+
+ match = {1., 1., 1., 1.};
+ return (size_t)-1;
+}
+
+
+template <typename T>
+size_t
+CPattern<T>::
+find( const SSignalRef<T>& signal,
+ ssize_t start,
+ int inc)
+{
+ if ( signal.samplerate != samplerate )
+ throw invalid_argument( "CPattern::find( SSignalRef&): not same samplerate");
+
+ return find( signal.signal,
+ start, inc);
+}
+
+template <typename T>
+size_t
+CPattern<T>::
+find( const valarray<T>& signal,
+ ssize_t start,
+ int inc)
+{
+ valarray<T> fenv_u, fenv_l;
+ envelope( {signal, samplerate}, SPatternPPack<T>::env_tightness,
+ 1./samplerate, &fenv_u, &fenv_l);
+
+ auto ftarget_freq =
+ exstrom::band_pass( signal, samplerate,
+ SPatternPPack<T>::bwf_ffrom,
+ SPatternPPack<T>::bwf_fupto,
+ SPatternPPack<T>::bwf_order, true);
+ auto fdzcdf =
+ dzcdf( SSignalRef<T> {signal, samplerate},
+ SPatternPPack<T>::dzcdf_step,
+ SPatternPPack<T>::dzcdf_sigma,
+ SPatternPPack<T>::dzcdf_smooth);
+
+ return find( fenv_u, fenv_l, ftarget_freq, fdzcdf,
+ start, inc);
+}
+
+
+// eof
diff --git a/src/sigproc/sigproc.cc b/src/sigproc/sigproc.cc
index a87e821..4266a3b 100644
--- a/src/sigproc/sigproc.cc
+++ b/src/sigproc/sigproc.cc
@@ -23,14 +23,11 @@ using namespace std;
template void sigproc::smooth( valarray<TFloat>&, size_t);
template void sigproc::normalize( valarray<TFloat>&);
template valarray<TFloat> sigproc::derivative( const valarray<TFloat>&);
-template size_t sigproc::envelope( const valarray<float>&, size_t, size_t, double, valarray<float>*, valarray<float>*, vector<size_t>*, vector<size_t>*);
-template size_t sigproc::envelope( const valarray<double>&, size_t, size_t, double, valarray<double>*, valarray<double>*, vector<size_t>*, vector<size_t>*);
-template valarray<TFloat> sigproc::dzcdf( const valarray<TFloat>&, size_t, float, float, size_t);
-template sigproc::CPattern<TFloat>::CPattern( const valarray<TFloat>&, size_t, size_t, size_t, const SPatternParamPack&, float, float, float);
-template size_t sigproc::CPattern<TFloat>::find( const valarray<TFloat>&, const valarray<TFloat>&, const valarray<TFloat>&, ssize_t, int);
-template size_t sigproc::CPattern<TFloat>::find( const valarray<TFloat>&, ssize_t, int);
+template size_t sigproc::envelope( const SSignalRef<float>&, size_t, double, valarray<float>*, valarray<float>*, vector<size_t>*, vector<size_t>*);
+template size_t sigproc::envelope( const SSignalRef<double>&, size_t, double, valarray<double>*, valarray<double>*, vector<size_t>*, vector<size_t>*);
+template valarray<TFloat> sigproc::dzcdf( const SSignalRef<TFloat>&, double, double, size_t);
template double sigproc::sig_diff( const valarray<TFloat>&, const valarray<TFloat>&, int);
-template double sigproc::phase_diff( const valarray<TFloat>&, const valarray<TFloat>&, size_t, size_t, size_t, float, float, unsigned, size_t);
+template double sigproc::phase_diff( const SSignalRef<TFloat>&, const SSignalRef<TFloat>&, size_t, size_t, double, double, unsigned, size_t);
valarray<float>
diff --git a/src/sigproc/sigproc.hh b/src/sigproc/sigproc.hh
index b875fb8..5f04a5c 100644
--- a/src/sigproc/sigproc.hh
+++ b/src/sigproc/sigproc.hh
@@ -5,7 +5,7 @@
* Author: Andrei Zavada <johnhommer at gmail.com>
* Initial version: 2011-01-26
*
- * Purpose: various standalone signal processing functions, class CPattern
+ * Purpose: various standalone signal processing functions
*
* License: GPL
*/
@@ -22,6 +22,7 @@
#include <samplerate.h>
#include "common/lang.hh"
+#include "common/alg.hh"
#include "exstrom.hh"
#if HAVE_CONFIG_H && !defined(VERSION)
@@ -32,6 +33,8 @@ using namespace std;
namespace sigproc {
+// simple functions operating irrespective of samplerate
+
template <typename T>
void
smooth( valarray<T>&, size_t side);
@@ -65,8 +68,6 @@ resample( const valarray<double>& signal,
int alg);
-
-
valarray<double>
interpolate_d( const vector<size_t>&,
unsigned, const valarray<double>&, float);
@@ -88,11 +89,22 @@ interpolate( const vector<size_t>& xi,
+
+// signal with samplerate
+
+template <typename T>
+struct SSignalRef {
+ const valarray<T>&
+ signal;
+ size_t samplerate;
+};
+
+
+
template <typename T>
size_t
-envelope( const valarray<T>& in,
+envelope( const SSignalRef<T>& in,
size_t dh, // tightness
- size_t samplerate,
double dt,
valarray<T>* env_lp = nullptr, // return interpolated
valarray<T>* env_up = nullptr, // return vector of points
@@ -100,79 +112,23 @@ envelope( const valarray<T>& in,
vector<size_t> *maxi_p = nullptr);
-
-
-
-template <typename T>
-int
-sign( const T& v)
-{
- return v >= 0. ? 1 : -1;
-}
-
-
-
template <typename T>
valarray<T>
-dzcdf( const valarray<T>& in,
- size_t samplerate,
- float dt,
- float sigma,
+dzcdf( const SSignalRef<T>& in,
+ double dt,
+ double sigma,
size_t smooth);
-
-template <typename T>
-struct SSignalRef {
- const valarray<T>&
- signal;
- unsigned
- samplerate;
-};
-
-
-
// cached signal property providers
template <typename T>
-struct SCachedDzcdf
- : public SSignalRef<T> {
- SCachedDzcdf (const valarray<T>& signal_, unsigned samplerate_)
- : SSignalRef<T> {signal_, samplerate_}
- {}
- SCachedDzcdf (const SCachedDzcdf&) = delete;
- // other ctors deleted implicitly due to a member of reference type
-
- const valarray<T>&
- operator()( float step_, float sigma_, unsigned smooth_)
- {
- if ( data.size() == 0 ||
- step != step_ || sigma != sigma_ || smooth != smooth_ )
- data = dzcdf<T>(
- SSignalRef<T>::signal, SSignalRef<T>::samplerate,
- step = step_, sigma = sigma_, smooth = smooth_);
- return data;
- }
- void drop()
- {
- data.resize(0);
- }
- private:
- float step,
- sigma;
- unsigned
- smooth;
- valarray<T>
- data;
-};
-
-template <typename T>
struct SCachedEnvelope
: public SSignalRef<T> {
- SCachedEnvelope (const valarray<T>& signal_, unsigned samplerate_)
- : SSignalRef<T> {signal_, samplerate_}
+ SCachedEnvelope (const SSignalRef<T>& signal_)
+ : SSignalRef<T> (signal_)
{}
SCachedEnvelope (const SCachedEnvelope&) = delete;
@@ -180,12 +136,14 @@ struct SCachedEnvelope
operator()( unsigned tightness_)
{
if ( lower.size() == 0 ||
- tightness != tightness_ )
- envelope( SSignalRef<T>::signal,
- tightness = tightness_, SSignalRef<T>::samplerate,
+ tightness != tightness_ ) {
+ envelope( (SSignalRef<T>)*this,
+ tightness = tightness_,
1./SSignalRef<T>::samplerate,
&lower,
&upper); // don't need anchor points, nor their count
+ mid = (upper + lower)/2;
+ }
return {lower, upper};
}
void drop()
@@ -194,7 +152,7 @@ struct SCachedEnvelope
lower.resize(0);
}
- float breadth( unsigned tightness_, size_t i)
+ T breadth( unsigned tightness_, size_t i)
{
(*this)( tightness_);
return upper[i] - lower[i];
@@ -205,24 +163,68 @@ struct SCachedEnvelope
return upper - lower;
}
+ T centre( unsigned tightness_, size_t i)
+ {
+ (*this)( tightness_);
+ return mid[i];
+ }
+ valarray<T> centre( unsigned tightness_)
+ {
+ (*this)( tightness_);
+ return mid;
+ }
+
private:
unsigned
tightness;
valarray<T>
upper,
+ mid,
lower;
};
template <typename T>
+struct SCachedDzcdf
+ : public SSignalRef<T> {
+ SCachedDzcdf (const SSignalRef<T>& signal_)
+ : SSignalRef<T> (signal_)
+ {}
+ SCachedDzcdf (const SCachedDzcdf&) = delete;
+ // other ctors deleted implicitly due to a member of reference type
+
+ const valarray<T>&
+ operator()( double step_, double sigma_, unsigned smooth_)
+ {
+ if ( data.size() == 0 ||
+ step != step_ || sigma != sigma_ || smooth != smooth_ )
+ data = dzcdf<T>(
+ (SSignalRef<T>)*this,
+ step = step_, sigma = sigma_, smooth = smooth_);
+ return data;
+ }
+ void drop()
+ {
+ data.resize(0);
+ }
+ private:
+ double step,
+ sigma;
+ unsigned
+ smooth;
+ valarray<T>
+ data;
+};
+
+template <typename T>
struct SCachedLowPassCourse
: public SSignalRef<T> {
- SCachedLowPassCourse (const valarray<T>& signal_, unsigned samplerate_)
- : SSignalRef<T> {signal_, samplerate_}
+ SCachedLowPassCourse (const SSignalRef<T>& signal_)
+ : SSignalRef<T> (signal_)
{}
SCachedLowPassCourse (const SCachedLowPassCourse&) = delete;
const valarray<T>&
- operator()( float fcutoff_, unsigned order_)
+ operator()( double fcutoff_, unsigned order_)
{
if ( data.size() == 0 ||
fcutoff != fcutoff_ || order != order_ )
@@ -237,7 +239,7 @@ struct SCachedLowPassCourse
}
private:
- float fcutoff;
+ double fcutoff;
unsigned
order;
valarray<TFloat>
@@ -247,13 +249,13 @@ struct SCachedLowPassCourse
template <typename T>
struct SCachedBandPassCourse
: public SSignalRef<T> {
- SCachedBandPassCourse (const valarray<T>& signal_, unsigned samplerate_)
- : SSignalRef<T> {signal_, samplerate_}
+ SCachedBandPassCourse (const SSignalRef<T>& signal_)
+ : SSignalRef<T> (signal_)
{}
SCachedBandPassCourse (const SCachedBandPassCourse&) = delete;
const valarray<T>&
- operator()( float ffrom_, float fupto_, unsigned order_)
+ operator()( double ffrom_, double fupto_, unsigned order_)
{
if ( data.size() == 0 ||
ffrom != ffrom_ || fupto != fupto_ || order != order_ )
@@ -268,7 +270,7 @@ struct SCachedBandPassCourse
}
private:
- float ffrom, fupto;
+ double ffrom, fupto;
unsigned
order;
valarray<TFloat>
@@ -279,102 +281,26 @@ struct SCachedBandPassCourse
-struct SPatternParamPack {
- int bwf_order;
- double bwf_ffrom,
- bwf_fupto;
- double dzcdf_step,
- dzcdf_sigma;
- int dzcdf_smooth,
- env_tightness;
- bool operator==( const SPatternParamPack& rv) const // cannot be defaulted!
- {
- return bwf_order == rv.bwf_order &&
- bwf_ffrom == rv.bwf_ffrom &&
- bwf_fupto == rv.bwf_fupto &&
- dzcdf_step == rv.dzcdf_step &&
- dzcdf_sigma == rv.dzcdf_sigma &&
- dzcdf_smooth == rv.dzcdf_smooth &&
- env_tightness == rv.env_tightness;
- }
-}; // keep fields in order, or edit ctor by initializer_list
-
-
-
-template <typename T>
-class CPattern {
- CPattern () = delete;
-
- public:
- // the complete pattern signature is made of:
- // (a) course of the mean (low-freq component);
- // (b) instantaneous frequency at fine intervals;
- // (c) signal breadth at given tightness.
-
- // data for individual constituents of the pattern:
- SPatternParamPack
- params;
-
- float a, b, c; // strictness
-
- // resulting
- float match_a,
- match_b,
- match_c;
-
- CPattern (const valarray<T>& pattern,
- size_t _context_before, size_t _context_after,
- size_t _samplerate,
- const SPatternParamPack&,
- float _a, float _b, float _c);
-
- size_t size_with_context() const
- {
- return course.size();
- }
- size_t size_essential() const
- {
- return size_with_context() - context_before - context_after;
- }
-
- size_t find( const valarray<T>& course,
- const valarray<T>& breadth,
- const valarray<T>& dzcdf,
- ssize_t start,
- int inc);
- size_t find( const valarray<T>& signal,
- ssize_t start,
- int inc);
-
- private:
- valarray<T>
- course,
- breadth,
- dzcd;
- size_t samplerate,
- context_before,
- context_after;
-};
-
-
template <typename T>
double
sig_diff( const valarray<T>& a, const valarray<T>& b, int d);
-
-
template <typename T>
double
-phase_diff( const valarray<T>& sig1,
- const valarray<T>& sig2,
- size_t samplerate,
+phase_diff( const SSignalRef<T>& sig1,
+ const SSignalRef<T>& sig2,
size_t sa, size_t sz,
- float fa, float fz,
+ double fa, double fz,
unsigned order,
size_t scope);
+
+
+
+
+
#include "sigproc.ii"
} // namespace sigproc
diff --git a/src/sigproc/sigproc.ii b/src/sigproc/sigproc.ii
index c937bd0..845ddb9 100644
--- a/src/sigproc/sigproc.ii
+++ b/src/sigproc/sigproc.ii
@@ -14,14 +14,11 @@ extern template void smooth( valarray<TFloat>&, size_t);
extern template void normalize( valarray<TFloat>&);
extern template valarray<TFloat> derivative( const valarray<TFloat>&);
// this one is used for both T = float and double
-extern template size_t envelope( const valarray<float>&, size_t, size_t, double, valarray<float>*, valarray<float>*, vector<size_t>*, vector<size_t>*);
-extern template size_t envelope( const valarray<double>&, size_t, size_t, double, valarray<double>*, valarray<double>*, vector<size_t>*, vector<size_t>*);
-extern template valarray<TFloat> dzcdf( const valarray<TFloat>&, size_t, float, float, size_t);
-extern template CPattern<TFloat>::CPattern( const valarray<TFloat>&, size_t, size_t, size_t, const SPatternParamPack&, float, float, float);
-extern template size_t CPattern<TFloat>::find( const valarray<TFloat>&, const valarray<TFloat>&, const valarray<TFloat>&, ssize_t, int);
-extern template size_t CPattern<TFloat>::find( const valarray<TFloat>&, ssize_t, int);
+extern template size_t envelope( const SSignalRef<float>&, size_t, double, valarray<float>*, valarray<float>*, vector<size_t>*, vector<size_t>*);
+extern template size_t envelope( const SSignalRef<double>&, size_t, double, valarray<double>*, valarray<double>*, vector<size_t>*, vector<size_t>*);
+extern template valarray<TFloat> dzcdf( const SSignalRef<TFloat>&, double, double, size_t);
extern template double sig_diff( const valarray<TFloat>&, const valarray<TFloat>&, int);
-extern template double phase_diff( const valarray<TFloat>&, const valarray<TFloat>&, size_t, size_t, size_t, float, float, unsigned, size_t);
+extern template double phase_diff( const SSignalRef<TFloat>&, const SSignalRef<TFloat>&, size_t, size_t, double, double, unsigned, size_t);
@@ -74,9 +71,8 @@ derivative( const valarray<T>& a)
template <typename T>
size_t
-envelope( const valarray<T>& in,
+envelope( const SSignalRef<T>& in,
size_t dh, // tightness
- size_t samplerate,
double dt,
valarray<T>* env_lp = nullptr, // return interpolated
valarray<T>* env_up = nullptr,
@@ -84,7 +80,7 @@ envelope( const valarray<T>& in,
vector<size_t> *maxi_p = nullptr)
{
size_t i, j,
- n_samples = in.size();
+ n_samples = in.signal.size();
vector<size_t>
mini,
@@ -96,19 +92,19 @@ envelope( const valarray<T>& in,
for ( i = dh; i < n_samples-dh; ++i ) {
for ( j = 1; j <= dh; ++j )
- if ( in[i-j] <= in[i] ) // [i] is not a local min
+ if ( in.signal[i-j] <= in.signal[i] ) // [i] is not a local min
goto inner_continue;
for ( j = 1; j <= dh; ++j )
- if ( in[i+j] <= in[i] ) // [i] is not
+ if ( in.signal[i+j] <= in.signal[i] ) // [i] is not
goto inner_continue;
mini.push_back( i);
continue;
inner_continue:
for ( j = 1; j <= dh; ++j )
- if ( in[i-j] >= in[i] ) // [i] is not a local max
+ if ( in.signal[i-j] >= in.signal[i] ) // [i] is not a local max
goto outer_continue;
for ( j = 1; j <= dh; ++j )
- if ( in[i+j] >= in[i] ) // [i] is not
+ if ( in.signal[i+j] >= in.signal[i] ) // [i] is not
goto outer_continue;
maxi.push_back( i);
outer_continue:
@@ -121,9 +117,9 @@ envelope( const valarray<T>& in,
if ( mini.size() > 5 && maxi.size() > 5 ) {
if ( env_lp )
- *env_lp = interpolate( mini, samplerate, in, dt);
+ *env_lp = interpolate( mini, in.samplerate, in.signal, dt);
if ( env_up )
- *env_up = interpolate( maxi, samplerate, in, dt);
+ *env_up = interpolate( maxi, in.samplerate, in.signal, dt);
if ( mini_p )
*mini_p = mini;
if ( maxi_p )
@@ -141,43 +137,38 @@ envelope( const valarray<T>& in,
template <typename T>
valarray<T>
-dzcdf( const valarray<T>& in,
- size_t samplerate,
- float dt,
- float sigma,
+dzcdf( const SSignalRef<T>& in,
+ double dt,
+ double sigma,
size_t smooth_side)
{
size_t i;
valarray<T>
- tmp (in),
- derivative (0., in.size());
+ tmp (in.signal),
+ deriv = derivative( in.signal);
smooth( tmp, smooth_side);
- // get derivative
- for ( i = 1; i < in.size(); ++i )
- derivative[i-1] = (tmp[i] - tmp[i-1]);
-
// collect zerocrossings
vector<size_t> izx;
- for ( i = 1; i < in.size(); ++i )
- if ( sign( derivative[i-1]) != sign( derivative[i]) )
+ for ( i = 1; i < in.signal.size(); ++i )
+ if ( agh::alg::sign( deriv[i-1]) != agh::alg::sign( deriv[i]) )
izx.push_back( i);
// prepare structures for interpolation
- size_t out_size = (float)in.size()/samplerate / dt;
+ size_t out_size = (double)in.signal.size()/in.samplerate / dt;
vector<size_t> xi (out_size);
- valarray<T> y (in.size());
+ valarray<T> y (in.signal.size());
// calculate the bloody zcdf
- float window = 4*dt; // half a second, good enough
- float t = 0., tdiff;
+ double window = 4*dt; // half a second, good enough
+ double t = 0., tdiff;
size_t I = 0, J;
for ( i = 0; i < out_size; ++i ) {
- xi[i] = i * dt * samplerate;
+ xi[i] = i * dt * in.samplerate;
for ( J = I; J > 0; --J ) {
- tdiff = (float)izx[J]/samplerate - t;
+ tdiff = (double)izx[J]/in.samplerate - t;
if ( tdiff > window )
continue;
if ( tdiff < -window )
@@ -185,7 +176,7 @@ dzcdf( const valarray<T>& in,
y[ xi[i] ] += exp( -gsl_pow_2(tdiff) / gsl_pow_2(sigma));
}
for ( ; J < izx.size(); ++J ) {
- tdiff = (float)izx[J]/samplerate - t;
+ tdiff = (double)izx[J]/in.samplerate - t;
if ( tdiff < -window )
continue;
if ( tdiff > window )
@@ -195,129 +186,7 @@ dzcdf( const valarray<T>& in,
t = i*dt;
I = J;
}
- return interpolate( xi, samplerate, y, 1./samplerate);
-}
-
-
-
-
-
-
-
-
-
-template <typename T>
-CPattern<T>::
-CPattern (const valarray<T>& pattern,
- size_t _context_before, size_t _context_after,
- size_t _samplerate,
- const SPatternParamPack& params_,
- float _a, float _b, float _c)
- : params (params_),
- a (_a), b (_b), c (_c),
- match_a (NAN), match_b (NAN), match_c (NAN),
- samplerate (_samplerate),
- context_before (_context_before), context_after (_context_after)
-{
- if ( context_before + context_after >= pattern.size() )
- throw invalid_argument ("pattern.size too small");
- course =
- exstrom::band_pass(
- pattern, samplerate,
- params.bwf_ffrom, params.bwf_fupto,
- params.bwf_order, true);
-
- valarray<T> env_u, env_l;
- envelope( pattern, params.env_tightness, samplerate,
- 1./samplerate,
- &env_l, &env_u);
- breadth.resize( env_u.size());
- breadth = env_u - env_l;
-
- dzcd = dzcdf( pattern, samplerate,
- params.dzcdf_step, params.dzcdf_sigma, params.dzcdf_smooth);
-}
-
-
-
-
-template <typename T>
-size_t
-CPattern<T>::
-find( const valarray<T>& fcourse,
- const valarray<T>& fbreadth,
- const valarray<T>& fdzcd,
- ssize_t start,
- int inc)
-{
- if ( inc == 0 || inc > (ssize_t)fcourse.size() ) {
- fprintf( stderr, "CSignalPattern::find(): bad search increment: %d\n", inc);
- return (size_t)-1;
- }
-
- T diff_course,
- diff_breadth,
- diff_dzcd;
-
- // printf( "course.size = %zu, fcourse.size = %zu, start = %zu\n",
- // course.size(), fcourse.size(), start);
- ssize_t iz = (inc > 0) ? fcourse.size() - size_with_context() : 0;
- size_t essential_part = size_essential();
- // bool looking_further = false;
- // T ax, bx, cx;
- for ( ssize_t i = start; (inc > 0) ? i < iz : i > iz; i += inc ) {
- diff_course = diff_breadth = diff_dzcd = 0.;
- for ( size_t j = 0; j < essential_part; ++j ) {
- diff_course += fdim( course [context_before + j], fcourse [i+j]);
- diff_breadth += fdim( breadth[context_before + j], fbreadth[i+j]);
- diff_dzcd += fdim( dzcd [context_before + j], fdzcd [i+j]);
- }
-
- diff_course /= essential_part;
- diff_breadth /= essential_part;
- diff_dzcd /= essential_part;
-
- // if ( i % 250 == 0 ) printf( "at %zu diff_course = %g,\tdiff_breadth = %g\t diff_dzcdf = %g\n", i, diff_course, diff_breadth, diff_dzcd);
- if ( diff_course < a && diff_breadth < b && diff_dzcd < c ) {
- // if ( !looking_further ) {
- // looking_further = true;
- match_a = diff_course, match_b = diff_breadth, match_c = diff_dzcd;
- return i;
- }
- }
-
- return (size_t)-1;
-}
-
-
-template <typename T>
-size_t
-CPattern<T>::
-find( const valarray<T>& signal,
- ssize_t start,
- int inc)
-{
- // low-pass signal being searched, too
- valarray<T> fcourse =
- exstrom::band_pass( signal, samplerate,
- params.bwf_ffrom, params.bwf_fupto,
- params.bwf_order, true);
-
- // prepare for comparison by other criteria:
- // signal envelope and breadth
- valarray<T> env_u, env_l;
- envelope( signal, params.env_tightness, samplerate,
- 1./samplerate, &env_u, &env_l);
- valarray<T> fbreadth (env_u.size());
- fbreadth = env_u - env_l;
-
- // dzcdf
- valarray<T> fdzcd =
- dzcdf( signal, samplerate,
- params.dzcdf_step, params.dzcdf_sigma, params.dzcdf_smooth);
-
- return find( fcourse, fbreadth, fdzcd,
- start, inc);
+ return interpolate( xi, in.samplerate, y, 1./in.samplerate);
}
@@ -341,20 +210,22 @@ sig_diff( const valarray<T>& a, const valarray<T>& b,
template <typename T>
double
-phase_diff( const valarray<T>& sig1,
- const valarray<T>& sig2,
- size_t samplerate,
+phase_diff( const SSignalRef<T>& sig1,
+ const SSignalRef<T>& sig2,
size_t sa, size_t sz,
- float fa, float fz,
+ double fa, double fz,
unsigned order,
size_t scope)
{
+ if ( sig1.samplerate != sig2.samplerate )
+ throw invalid_argument ("sigproc::phase_diff(): sig1.samplerate != sig2.samplerate");
if ( order == 0 )
- throw invalid_argument ("NExstrom::phase_diff(): order == 0");
+ throw invalid_argument ("sigproc::phase_diff(): order == 0");
+
// bandpass sig1 and sig2
valarray<T>
- sig1p = exstrom::band_pass( valarray<T> (&sig1[sa], sz - sa), samplerate, fa, fz, order, true),
- sig2p = exstrom::band_pass( valarray<T> (&sig2[sa], sz - sa), samplerate, fa, fz, order, true);
+ sig1p = exstrom::band_pass( valarray<T> (&sig1.signal[sa], sz - sa), sig1.samplerate, fa, fz, order, true),
+ sig2p = exstrom::band_pass( valarray<T> (&sig2.signal[sa], sz - sa), sig2.samplerate, fa, fz, order, true);
// slide one against the other a little
double diff = INFINITY, old_diff, diff_min = INFINITY;
@@ -374,7 +245,7 @@ phase_diff( const valarray<T>& sig1,
diff_min = diff, dist_min = dist;
} while ( (dist++) < (int)scope && old_diff > diff );
- return (double)dist_min / samplerate;
+ return (double)dist_min / sig1.samplerate;
}
diff --git a/src/ui/sf/sf-channel.cc b/src/ui/sf/sf-channel.cc
index cf928ba..51fad68 100644
--- a/src/ui/sf/sf-channel.cc
+++ b/src/ui/sf/sf-channel.cc
@@ -35,10 +35,10 @@ SChannel( agh::CRecording& r,
annotations (r.F().annotations(name)),
artifacts (r.F().artifacts(name)),
_p (parent),
- signal_lowpass (signal_filtered, samplerate()),
- signal_bandpass (signal_filtered, samplerate()),
- signal_envelope (signal_filtered, samplerate()),
- signal_dzcdf (signal_filtered, samplerate()),
+ signal_lowpass ({signal_filtered, samplerate()}),
+ signal_bandpass ({signal_filtered, samplerate()}),
+ signal_envelope ({signal_filtered, samplerate()}),
+ signal_dzcdf ({signal_filtered, samplerate()}),
zeroy (y0),
// let them be read or recalculated
signal_display_scale( NAN),
@@ -138,8 +138,8 @@ SChannel( agh::CRecording& r,
} else if ( type == sigfile::SChannel::TType::emg ) {
valarray<TFloat> env_u, env_l;
- sigproc::envelope( signal_original,
- 5, samplerate(), 1.,
+ sigproc::envelope( {signal_original, samplerate()},
+ 5, 1.,
&env_l, &env_u);
emg_profile.resize( env_l.size());
emg_profile = env_u - env_l;
diff --git a/src/ui/sf/sf.hh b/src/ui/sf/sf.hh
index 0548cfd..ad55f71 100644
--- a/src/ui/sf/sf.hh
+++ b/src/ui/sf/sf.hh
@@ -24,6 +24,7 @@
#include "common/config-validate.hh"
#include "sigproc/winfun.hh"
#include "sigproc/sigproc.hh"
+#include "sigproc/patterns.hh"
#include "metrics/phasic-events.hh"
#include "expdesign/primaries.hh"
#include "ica/ica.hh"
@@ -501,7 +502,7 @@ class SScoringFacility
DELETE_DEFAULT_METHODS (SFindDialog);
// own copies of parent's same
- sigproc::SPatternParamPack
+ sigproc::SPatternPPack<TFloat>
Pp,
Pp2;
--
Sleep experiment manager
More information about the debian-med-commit
mailing list