# HG changeset patch
# User Chris Cannam
# Date 1439228661 -3600
# Node ID d0880801415d23079214d3f08c5a8caf77502fb4
# Parent  f198bfdf036698c7dc2d0f9307760015dd62c0c8
Filter code, delays

diff -r f198bfdf0366 -r d0880801415d calc_delays.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/calc_delays.m	Mon Aug 10 18:44:21 2015 +0100
@@ -0,0 +1,14 @@
+
+load MIDI_FB_ellip_pitch_60_96_22050_Q25.mat;
+
+dirac = zeros(50000, 1);
+dirac(1) = 1.0;
+
+delays = zeros(120, 1);
+
+for n = 21:120
+    f = filter(h(n).b, h(n).a, dirac);
+    [~,pos] = max(f(1:10000));
+    [~,neg] = max(-f(1:10000));
+    delays(n) = 1 + (pos + neg) / 2;
+end
diff -r f198bfdf0366 -r d0880801415d delays.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/delays.h	Mon Aug 10 18:44:21 2015 +0100
@@ -0,0 +1,123 @@
+static const int filter_delay[120] = {
+0,
+0,
+0,
+0,
+0,
+0,
+0,
+0,
+0,
+0,
+0,
+0,
+0,
+0,
+0,
+0,
+0,
+0,
+0,
+0,
+459,
+434,
+395,
+387,
+709,
+681,
+654,
+618,
+583,
+541,
+511,
+482,
+455,
+430,
+399,
+383,
+362,
+341,
+323,
+299,
+293,
+272,
+257,
+242,
+237,
+212,
+211,
+193,
+185,
+172,
+157,
+156,
+137,
+142,
+127,
+122,
+113,
+107,
+107,
+487,
+467,
+434,
+403,
+380,
+353,
+333,
+326,
+302,
+291,
+274,
+246,
+240,
+219,
+218,
+195,
+191,
+184,
+168,
+167,
+152,
+144,
+138,
+133,
+121,
+124,
+108,
+99,
+97,
+88,
+92,
+82,
+77,
+74,
+83,
+80,
+305,
+298,
+263,
+244,
+234,
+225,
+213,
+201,
+190,
+179,
+152,
+163,
+156,
+143,
+135,
+132,
+118,
+131,
+111,
+100,
+93,
+91,
+86,
+80,
+79
+};
+
diff -r f198bfdf0366 -r d0880801415d src/Filter.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Filter.cpp	Mon Aug 10 18:44:21 2015 +0100
@@ -0,0 +1,66 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#include "Filter.h"
+
+#include <bqvec/VectorOps.h>
+#include <bqvec/Allocators.h>
+
+#include <stdexcept>
+
+using namespace std;
+using namespace breakfastquay;
+
+Filter::Filter(Parameters params)
+{
+    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_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_inbuf = allocate_and_zero<double>(m_sz);
+    m_outbuf = allocate_and_zero<double>(m_order);
+}
+
+Filter::~Filter()
+{
+    deallocate(m_a);
+    deallocate(m_b);
+    deallocate(m_inbuf);
+    deallocate(m_outbuf);
+}
+
+void
+Filter::reset()
+{
+    v_zero(m_inbuf, m_sz);
+    v_zero(m_outbuf, m_order);
+}
+
+void
+Filter::process(const double *const BQ_R__ in,
+		double *const BQ_R__ out,
+		const int n)
+{
+    for (int s = 0; s < n; ++s) {
+
+	v_move(m_inbuf + 1, m_inbuf, m_order);
+	m_inbuf[0] = in[s];
+
+	double b_sum = v_multiply_and_sum(m_b, m_inbuf, m_sz);
+	double a_sum = v_multiply_and_sum(m_a + 1, m_outbuf, m_order);
+
+	double outval = b_sum - a_sum;
+
+	v_move(m_outbuf + 1, m_outbuf, m_order - 1);
+	m_outbuf[0] = outval;
+	out[s] = outval;
+    }
+}
+
diff -r f198bfdf0366 -r d0880801415d src/Filter.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/Filter.h	Mon Aug 10 18:44:21 2015 +0100
@@ -0,0 +1,35 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+#ifndef FILTER_H
+#define FILTER_H
+
+#include <bqvec/Restrict.h>
+
+#include <vector>
+
+class Filter
+{
+public:
+    struct Parameters {
+        std::vector<double> a;
+        std::vector<double> b;
+    };
+
+    Filter(Parameters params);
+    ~Filter();
+
+    void reset();
+    void process(const double *const BQ_R__ in,
+                 double *const BQ_R__ out,
+                 const int n);
+
+private:
+    int m_order;
+    int m_sz;
+    double *m_a;
+    double *m_b;
+    double *m_inbuf;
+    double *m_outbuf;
+};
+    
+#endif
diff -r f198bfdf0366 -r d0880801415d src/test-filter.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/test-filter.cpp	Mon Aug 10 18:44:21 2015 +0100
@@ -0,0 +1,36 @@
+
+#include "Filter.h"
+
+#include <iostream>
+#include <cmath>
+
+using namespace std;
+
+int main(int argc, char **argv)
+{
+    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 };
+    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 };
+
+    bool good = true;
+    double thresh = 1e-14;
+    
+    for (int i = 0; i < 8; ++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;
+}
+