Mercurial > hg > qm-dsp
changeset 506:285f18c0992a
Tests and fixes for FiltFilt: Fix overrun; reset filter between forwards and backwards processes
| author | Chris Cannam <cannam@all-day-breakfast.com> |
|---|---|
| date | Wed, 05 Jun 2019 15:50:38 +0100 |
| parents | 930b5b0f707d |
| children | d7b9691817a3 |
| files | dsp/signalconditioning/FiltFilt.cpp dsp/signalconditioning/FiltFilt.h tests/Makefile tests/TestFiltFilt.cpp |
| diffstat | 4 files changed, 127 insertions(+), 17 deletions(-) [+] |
line wrap: on
line diff
--- a/dsp/signalconditioning/FiltFilt.cpp Wed Jun 05 12:55:15 2019 +0100 +++ b/dsp/signalconditioning/FiltFilt.cpp Wed Jun 05 15:50:38 2019 +0100 @@ -15,10 +15,6 @@ #include "FiltFilt.h" -////////////////////////////////////////////////////////////////////// -// Construction/Destruction -////////////////////////////////////////////////////////////////////// - FiltFilt::FiltFilt(Filter::Parameters parameters) : m_filter(parameters) { @@ -29,7 +25,9 @@ { } -void FiltFilt::process(double *src, double *dst, int length) +void FiltFilt::process(const double *const QM_R__ src, + double *const QM_R__ dst, + const int length) { int i; @@ -53,19 +51,25 @@ int index = 0; for (i = nFact; i > 0; i--) { - filtScratchIn[ index++ ] = sample0 - src[ i ]; + if (i < length) { + filtScratchIn[index] = sample0 - src[ i ]; + } + ++index; } index = 0; for (i = 0; i < nFact; i++) { - filtScratchIn[ (nExt - nFact) + index++ ] = - sampleN - src[ (length - 2) - i ]; + if (i < length) { + filtScratchIn[(nExt - nFact) + index] = + sampleN - src[ (length - 2) - i ]; + } + ++index; } index = 0; for (i = 0; i < length; i++) { filtScratchIn[ i + nFact ] = src[ i ]; } - + //////////////////////////////// // Do 0Ph filtering m_filter.process(filtScratchIn, filtScratchOut, nExt); @@ -75,9 +79,12 @@ filtScratchIn[ i ] = filtScratchOut[ nExt - i - 1]; } + // clear filter state + m_filter.reset(); + // do FILTER again m_filter.process(filtScratchIn, filtScratchOut, nExt); - + // reverse the series back for (i = 0; i < nExt; i++) { filtScratchIn[ i ] = filtScratchOut[ nExt - i - 1 ]; @@ -95,7 +102,3 @@ delete [] filtScratchOut; } -void FiltFilt::reset() -{ - -}
--- a/dsp/signalconditioning/FiltFilt.h Wed Jun 05 12:55:15 2019 +0100 +++ b/dsp/signalconditioning/FiltFilt.h Wed Jun 05 15:50:38 2019 +0100 @@ -28,8 +28,9 @@ FiltFilt(Filter::Parameters); virtual ~FiltFilt(); - void reset(); - void process(double* src, double* dst, int length); + void process(const double *const QM_R__ src, + double *const QM_R__ dst, + const int length); private: Filter m_filter;
--- a/tests/Makefile Wed Jun 05 12:55:15 2019 +0100 +++ b/tests/Makefile Wed Jun 05 15:50:38 2019 +0100 @@ -5,7 +5,7 @@ LDFLAGS := $(LDFLAGS) -lboost_unit_test_framework -lpthread LIBS := ../libqm-dsp.a -TESTS := test-mathutilities test-window test-filter test-fft test-dct test-pvoc test-resampler test-medianfilter test-getkeymode test-chromagram +TESTS := test-mathutilities test-window test-filter test-filtfilt test-fft test-dct test-pvoc test-resampler test-medianfilter test-getkeymode test-chromagram VG := valgrind -q @@ -25,6 +25,9 @@ test-filter: TestFilter.o $(LIBS) $(CXX) -o $@ $^ $(LDFLAGS) +test-filtfilt: TestFiltFilt.o $(LIBS) + $(CXX) -o $@ $^ $(LDFLAGS) + test-fft: TestFFT.o $(LIBS) $(CXX) -o $@ $^ $(LDFLAGS) @@ -47,6 +50,7 @@ TestMedianFilter.o: $(LIBS) TestWindow.o: $(LIBS) TestFilter.o: $(LIBS) +TestFiltFilt.o: $(LIBS) TestFFT.o: $(LIBS) TestDCT.o: $(LIBS) TestPhaseVocoder.o: $(LIBS)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/TestFiltFilt.cpp Wed Jun 05 15:50:38 2019 +0100 @@ -0,0 +1,102 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ + +#include "dsp/signalconditioning/FiltFilt.h" + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MAIN + +#include <boost/test/unit_test.hpp> + +#include <stdexcept> + +using namespace std; + +BOOST_AUTO_TEST_SUITE(TestFiltFilt) + +// Our FiltFilt is, apparently intentionally (?) a bit different from +// the equivalents in MATLAB and Octave (which themselves also differ) +// in how it handles transients at either end. While the three should +// all be "roughly" the same, all produce slightly different values. +// So the "expected" values here are, unfortunately, simply produced +// by the function itself, in a process involving tracing the values +// passed in to each of the filter stages within it and compared their +// outputs against MATLAB. + +static vector<double> iir_a { 1,5.75501989315662,16.326056867468,28.779190797823,34.2874379215653,28.137815126537,15.6064643257793,5.37874515231553,0.913800050254382,0.0,0.0 }; +static vector<double> iir_b { 0.0031954608137085,0.0180937089815597,0.0508407778575426,0.0895040074158415,0.107385387168148,0.0895040074158415,0.0508407778575426,0.0180937089815597,0.0031954608137085,0.0,0.0 }; + +BOOST_AUTO_TEST_CASE(long_zeroes) +{ + vector<double> a(iir_a); + vector<double> b(iir_b); + + vector<double> zeroes(1000, 0.0); + vector<double> expected(zeroes); + + FiltFilt f({ a, b }); + + int n = expected.size(); + vector<double> out(n, 0.0); + + f.process(zeroes.data(), out.data(), n); + + double thresh = 1e-12; + + for (int i = 0; i < n; ++i) { + BOOST_CHECK_SMALL(out[i] - expected[i], thresh); + } +} + +BOOST_AUTO_TEST_CASE(empty) +{ + vector<double> a(iir_a); + vector<double> b(iir_b); + vector<double> in; + vector<double> out; + vector<double> expected; + + FiltFilt f({ a, b }); + + f.process(in.data(), out.data(), 0); +} + +BOOST_AUTO_TEST_CASE(single_value) +{ + vector<double> a(iir_a); + vector<double> b(iir_b); + vector<double> in({ 1.0 }); + vector<double> out({ 0.0 }); + vector<double> expected({ -0.000236642942007 }); + + FiltFilt f({ a, b }); + + f.process(in.data(), out.data(), 1); + + double thresh = 1e-12; + + BOOST_CHECK_SMALL(out[0] - expected[0], thresh); +} + +BOOST_AUTO_TEST_CASE(shortish) +{ + vector<double> a(iir_a); + vector<double> b(iir_b); + + vector<double> in { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 }; + vector<double> expected { 0.0180395628158, -0.0317213769922, 0.0276928584707, -0.00944844359572, -0.0119394068579, 0.0252763627734, -0.0239770024451, 0.010560830803, 0.00726290602379, -0.0191915890394, 0.0202451156989, -0.0103434265023, -0.00341952789985, 0.0142823904341, -0.0160627659442, 0.00984876782803 }; + + FiltFilt f({ a, b }); + + int n = expected.size(); + vector<double> out(n, 0.0); + + f.process(in.data(), out.data(), n); + + double thresh = 1e-12; + + for (int i = 0; i < n; ++i) { + BOOST_CHECK_SMALL(out[i] - expected[i], thresh); + } +} + +BOOST_AUTO_TEST_SUITE_END()