[med-svn] [SCM] aghermann branch, master, updated. 551e213a23b59b71cba6a9c3a282d1b60e21b854

[Andrei Zavada]

The following commit has been merged in the master branch:
commit 081c27332383f295e2518ef5b96055a8b6f087f3
Author: Andrei Zavada <johnhommer at gmail.com>
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