The following commit has been merged in the master branch:
commit 081c27332383f295e2518ef5b96055a8b6f087f3
Author: Andrei Zavada <[email protected]>
Date: Mon Apr 22 01:17:53 2013 +0300
refactor ext-filters, fix one uninitialised, but useful, read
diff --git a/src/sigproc/ext-filters.hh b/src/sigproc/ext-filters.hh
index 836d55c..85c598f 100644
--- a/src/sigproc/ext-filters.hh
+++ b/src/sigproc/ext-filters.hh
@@ -130,26 +130,24 @@ class CFilterSE
T ts = 1.0 / CFilterIIR<T>::samplerate,
fprewarp, r, s, t;
- fprewarp = tan(f0 * M_PI * ts) / (M_PI * ts);
- r = gsl_pow_2(2. * M_PI * fprewarp * ts);
+ fprewarp = tan( f0 * M_PI * ts) / (M_PI * ts);
+ r = gsl_pow_2( 2. * M_PI * fprewarp * ts);
// From November 1992 prewarping applied because of
Arends results !
// r:=sqr(2.0*pi*f0*Ts); No
prewarping
s = 2. * M_PI * bandwidth * ts * 2.;
t = 4. + r + s;
- CFilterIIR<T>::poles =
- { (T)1.,
- (T)((8.0 - 2.0 * r) / t),
- (T)((-4.0 + s - r) / t) };
+ CFilterIIR<T>::poles[0] = 1.;
+ CFilterIIR<T>::poles[1] = (8.0 - 2.0 * r) / t;
+ CFilterIIR<T>::poles[2] = (-4.0 + s - r) / t;
fprewarp = tan(fc * M_PI * ts) / (M_PI * ts);
r = 2.0 / (2. * M_PI * fprewarp);
s = CFilterIIR<T>::gain * 2. * M_PI * bandwidth * 2.;
- CFilterIIR<T>::zeros =
- { (T)(s * (r + ts) / t),
- (T)(s * (-2.0 * r) / t),
- (T)(s * (r - ts) / t) };
- }
+ CFilterIIR<T>::zeros[0] = s * (r + ts) / t;
+ CFilterIIR<T>::zeros[1] = s * (-2.0 * r) / t;
+ CFilterIIR<T>::zeros[2] = s * (r - ts) / t;
+ }
private:
T f0,
@@ -169,7 +167,7 @@ class CFilterDUE
: CFilterIIR<T> (samplerate_, direction_, gain_, back_polate_),
minus_3db_frequency (minus_3db_frequency_)
{
- CFilterIIR<T>::zeros.resize(2);
CFilterIIR<T>::filter_state_z.resize(2);
+ CFilterIIR<T>::zeros.resize(3);
CFilterIIR<T>::filter_state_z.resize(3);
CFilterIIR<T>::poles.resize(1);
CFilterIIR<T>::filter_state_p.resize(2); // NrPoles+1 !!!!!
calculate_iir_coefficients();
}
@@ -182,12 +180,25 @@ class CFilterDUE
fprewarp = tan( M_PI * minus_3db_frequency *
ts) / (M_PI * ts),
r = 1. / (2. * M_PI * fprewarp),
s = ts / 2.;
- CFilterIIR<T>::zeros = {
- (CFilterIIR<T>::gain * (s + r)),
- (CFilterIIR<T>::gain * (s - r)),
- 1.};
- }
+ /// this is what m.roussen has in
Library/Filters/DUEFilter.cs:
+ // FZeros[0] = Gain * (s + r);
+ // FZeros[1] = Gain * (s - r);
+ // FPoles[0] = 1.0;
+ /// note the last assignment is to FPoles -- which
isn't used anyway at index 0
+ CFilterIIR<T>::zeros[0] = CFilterIIR<T>::gain * (s + r);
+ CFilterIIR<T>::zeros[1] = CFilterIIR<T>::gain * (s - r);
+ CFilterIIR<T>::zeros[2] = 1.;
+
+ // so I got zeros[2] assigned by transcription mistake
(should be to poles[0]) -- and
+ // it worked! *Except* that zeros, through this
assignment, has now grown to 3, which
+ // causes an invalid read in apply(), and went
miraculously undetected until noticed
+ // (thank valgrind) at 0.9_rc stage
+ CFilterIIR<T>::poles[0] = 1.;
+ // Still, FPoles[0] = 1.0 is here for whatever it is
intended to do
+
+ // May your life be forever geweldig, mate!
+ }
private:
T minus_3db_frequency;
diff --git a/src/sigproc/ext-filters.ii b/src/sigproc/ext-filters.ii
index 7ba1863..581e251 100644
--- a/src/sigproc/ext-filters.ii
+++ b/src/sigproc/ext-filters.ii
@@ -16,7 +16,7 @@ extern template valarray<double> CFilterIIR<double>::apply(
const valarray<doubl
template <typename T>
valarray<T>
-sigproc::CFilterIIR<T>::
+CFilterIIR<T>::
apply( const valarray<T>& in, bool use_first_sample_to_reset)
{
if ( unlikely (poles.size() == 0) )
@@ -49,29 +49,25 @@ apply( const valarray<T>& in, bool
use_first_sample_to_reset)
use_first_sample_to_reset = false;
}
// Compute new output sample
- double r = 0.;
+ T R = 0.;
// Add past output-values
size_t j;
for ( j = 1; j < poles.size(); ++j )
- r += poles[j] * filter_state_p[j];
+ R += poles[j] * filter_state_p[j];
// Not anticipate = do not include current input sample in
output value
- if ( not anticipate)
- out[i] = r;
- // Add past input-values
- for ( j = 0; j < zeros.size(); ++j )
- r += zeros[j] * filter_state_z[j];
// Anticipate = include current input sample in output value
- if ( anticipate )
- out[i] = r;
+ if ( anticipate) // Add past input-values
+ for ( j = 0; j < zeros.size(); ++j )
+ R += zeros[j] * filter_state_z[j];
// Do backpolation (FilterStateP[1] = Last out-sample)
- out[i] = back_polate * filter_state_p[1] + (1.0 - back_polate)
* out[i];
+ out[i] = back_polate * filter_state_p[1] + (1.0 - back_polate)
* R;
// Scale result
// TODO: Check if removing extra checks was ok
// Update filter state
- for ( j = poles.size()-1; j >= 2; --j )
+ for ( j = filter_state_p.size()-1; j >= 2; --j )
filter_state_p[j] = filter_state_p[j-1];
- filter_state_p[1] = r;
- for ( j = zeros.size()-1; j >= 1; --j )
+ filter_state_p[1] = R;
+ for ( j = filter_state_z.size()-1; j >= 1; --j )
filter_state_z[j] = filter_state_z[j-1];
}
--
Sleep experiment manager
_______________________________________________
debian-med-commit mailing list
[email protected]
http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
