Mercurial > hg > tipic

changeset 34:156d560d3470

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Support FIR filters in Filter
author Chris Cannam
date Wed, 30 Sep 2015 15:27:04 +0100
parents e5d5a7098a32
children 474f45cfd3fd
files src/Filter.cpp src/Filter.h src/test-filter.cpp
diffstat 3 files changed, 82 insertions(+), 25 deletions(-) [+]
line wrap: on
line diff
--- a/src/Filter.cpp	Wed Sep 30 15:08:08 2015 +0100
+++ b/src/Filter.cpp	Wed Sep 30 15:27:04 2015 +0100
@@ -12,31 +12,45 @@
 
 Filter::Filter(Parameters params)
 {
-    if (params.a.size() != params.b.size()) {
-	throw logic_error("Inconsistent numbers of filter coefficients");
+    if (params.a.empty()) {
+        m_fir = true;
+        if (params.b.empty()) {
+            throw logic_error("Filter must have at least one pair of coefficients");
+        }
+    } else {
+        m_fir = false;
+        if (params.a.size() != params.b.size()) {
+            throw logic_error("Inconsistent numbers of filter coefficients");
+        }
     }
-    if (params.a.size() < 1) {
-	throw logic_error("Filter must have at least one pair of coefficients");
-    }
-    m_sz = int(params.a.size());
+    
+    m_sz = int(params.b.size());
     m_order = m_sz - 1;
-    m_a = allocate<double>(m_sz);
-    m_b = allocate<double>(m_sz);
-    v_copy(m_a, &params.a[0], m_sz);
-    v_copy(m_b, &params.b[0], m_sz);
+
     m_offmax = 20;
     m_offa = m_offmax;
     m_offb = m_offmax;
-    m_bufa = allocate_and_zero<double>(m_order + m_offmax);
+
+    if (m_fir) {
+        m_a = 0;
+        m_bufa = 0;
+    } else {
+        m_a = allocate<double>(m_sz);
+        m_bufa = allocate_and_zero<double>(m_order + m_offmax);
+        v_copy(m_a, &params.a[0], m_sz);
+    }
+    
+    m_b = allocate<double>(m_sz);
+    v_copy(m_b, &params.b[0], m_sz);
     m_bufb = allocate_and_zero<double>(m_sz + m_offmax);
 }
 
 Filter::~Filter()
 {
     deallocate(m_a);
+    deallocate(m_bufa);
     deallocate(m_b);
     deallocate(m_bufb);
-    deallocate(m_bufa);
 }
 
 void
@@ -45,7 +59,7 @@
     m_offb = m_offmax;
     m_offa = m_offmax;
     v_zero(m_bufb, m_sz + m_offmax);
-    v_zero(m_bufa, m_order + m_offmax);
+    if (!m_fir) v_zero(m_bufa, m_order + m_offmax);
 }
 
 void
@@ -63,16 +77,25 @@
         m_bufb[m_offb] = in[s];
 
 	double b_sum = v_multiply_and_sum(m_b, m_bufb + m_offb, m_sz);
-	double a_sum = v_multiply_and_sum(m_a + 1, m_bufa + m_offa, m_order);
 
-	double outval = b_sum - a_sum;
+        double outval;
 
-        if (m_offa > 0) --m_offa;
-        else {
-            v_move(m_bufa + m_offmax + 1, m_bufa, m_order - 1);
-            m_offa = m_offmax;
+        if (m_fir) {
+
+            outval = b_sum;
+
+        } else {
+
+            double a_sum = v_multiply_and_sum(m_a + 1, m_bufa + m_offa, m_order);
+            outval = b_sum - a_sum;
+
+            if (m_offa > 0) --m_offa;
+            else {
+                v_move(m_bufa + m_offmax + 1, m_bufa, m_order - 1);
+                m_offa = m_offmax;
+            }
+            m_bufa[m_offa] = outval;
         }
-        m_bufa[m_offa] = outval;
         
 	out[s] = outval;
     }
--- a/src/Filter.h	Wed Sep 30 15:08:08 2015 +0100
+++ b/src/Filter.h	Wed Sep 30 15:27:04 2015 +0100
@@ -15,7 +15,14 @@
         std::vector<double> b;
     };
 
+    /**
+     * Construct an IIR filter with numerators b and denominators
+     * a. To make an FIR filter instead (i.e. with denominator = 1),
+     * leave the vector a in the param struct empty. Otherwise, the
+     * vectors a and b must have the same number of values.
+     */
     Filter(Parameters params);
+    
     ~Filter();
 
     void reset();
@@ -39,6 +46,7 @@
     int m_offa;
     int m_offb;
     int m_offmax;
+    bool m_fir;
 
     Filter(const Filter &); // not supplied
     Filter &operator=(const Filter &); // not supplied
--- a/src/test-filter.cpp	Wed Sep 30 15:08:08 2015 +0100
+++ b/src/test-filter.cpp	Wed Sep 30 15:27:04 2015 +0100
@@ -8,22 +8,26 @@
 
 int main(int argc, char **argv)
 {
+    // IIR
+    
     vector<double> a { 1,5.75501989315662,16.326056867468,28.779190797823,34.2874379215653,28.137815126537,15.6064643257793,5.37874515231553,0.913800050254382,0.0,0.0 };
     vector<double> b { 0.0031954608137085,0.0180937089815597,0.0508407778575426,0.0895040074158415,0.107385387168148,0.0895040074158415,0.0508407778575426,0.0180937089815597,0.0031954608137085,0.0,0.0 };
 
     Filter f({ a, b });
     
     vector<double> in { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 };
-    vector<double> out(8, 0.0);
-
-    f.process(in.data(), out.data(), 8);
 
     vector<double> expected { 0.003195460813709, 0.006094690058282, 0.009370240771381, 0.012857578361690, 0.015328760300750, 0.019107809614909, 0.022257958968869, 0.024598034053011, 0.029106103380941, 0.031152166476509, 0.034424013713795, 0.038775350541015, 0.039924063374886, 0.044846280036012, 0.047614917256999, 0.049338485830505 };
+
+    int n = expected.size();
+    vector<double> out(n, 0.0);
+
+    f.process(in.data(), out.data(), n);
     
     bool good = true;
-    double thresh = 1e-14;
+    double thresh = 1e-12;
     
-    for (int i = 0; i < 8; ++i) {
+    for (int i = 0; i < n; ++i) {
 	if (fabs(out[i] - expected[i]) > thresh) {
 	    cerr << "ERROR: out[" << i << "] (" << out[i]
 		 << ") differs from expected[" << i << "] (" << expected[i]
@@ -32,6 +36,28 @@
 	}
     }
 
+    // FIR
+
+    b = { -1.5511e-18,-0.022664,1.047e-17,0.27398,0.49737,0.27398,1.047e-17,-0.022664,-1.5511e-18 };
+    Filter ff({ {}, b });
+
+    expected = { -1.5511e-18,-0.022664,-0.045328,0.20599,0.95467,1.9773,3,4,5,6,7,8,9,10,11,12 };
+
+    n = expected.size();
+    
+    ff.process(in.data(), out.data(), n);
+
+    thresh = 1e-4;
+    
+    for (int i = 0; i < n; ++i) {
+	if (fabs(out[i] - expected[i]) > thresh) {
+	    cerr << "ERROR: out[" << i << "] (" << out[i]
+		 << ") differs from expected[" << i << "] (" << expected[i]
+		 << ") by " << out[i] - expected[i] << endl;
+	    good = false;
+	}
+    }
+    
     if (good) {
 	cerr << "Success" << endl;
 	return 0;