diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..cd71fc8 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +/build/ +/.vscode/ +/build-*/ \ No newline at end of file diff --git a/Biquad.cpp b/Biquad.cpp index ba22b07..3957a80 100644 --- a/Biquad.cpp +++ b/Biquad.cpp @@ -36,6 +36,12 @@ Biquad::Biquad(){ } + +Biquad::Biquad(double b0, double b1, double b2, double a1, double a2): + b0(b0),b1(b1),b2(b2),a1(a1),a2(a2) +{ +} + Biquad::~Biquad(){ } diff --git a/Biquad.h b/Biquad.h index c4f63a1..f1146ed 100644 --- a/Biquad.h +++ b/Biquad.h @@ -27,6 +27,8 @@ */ +#ifndef BIQUAD_H +#define BIQUAD_H #include #include @@ -36,10 +38,13 @@ using namespace std; // A biquad filter expression: // y[n] = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] + a1 * y[n-1] + a2 * y[n-2]; +// DF2T: Direct Form II Transposed, It is one of the standard structures used to implement digital filters, including Butterworth filters. + class Biquad { public: Biquad(); + Biquad(double b0, double b1, double b2, double a1, double a2); ~Biquad(); @@ -93,3 +98,5 @@ class BiquadChain { void allocate(int count); }; + +#endif \ No newline at end of file diff --git a/Butterworth.cpp b/Butterworth.cpp index 9a7a672..180c2b9 100644 --- a/Butterworth.cpp +++ b/Butterworth.cpp @@ -28,7 +28,7 @@ */ -#import +#include #include #include "Butterworth.h" diff --git a/Butterworth.h b/Butterworth.h index 9debe3b..cddd486 100644 --- a/Butterworth.h +++ b/Butterworth.h @@ -27,6 +27,8 @@ */ +#ifndef BUTTERWORTH_H +#define BUTTERWORTH_H #include "Biquad.h" @@ -116,3 +118,4 @@ class Butterworth { double * ba; }; +#endif diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..dec2278 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,36 @@ +# 交叉编译为 arm 架构程序。注意: 本块必须放在最前面，否则不生效 +if(ARM64) + message(STATUS "Using ARM64 Tool Chain!") + set(CMAKE_TOOLCHAIN_FILE ${CMAKE_SOURCE_DIR}/toolchain-arm64.cmake) + # 总是 Release 版本 + set(CMAKE_BUILD_TYPE Release) + message(STATUS "Compiler:${CMAKE_CXX_COMPILER}") +else() + message(STATUS "Using X64 Tool Chain!") + # 指定编译器选项 + # 为了 gdb 调试，不优化、带上debug符号 + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O0 -g -rdynamic") +endif() + +cmake_minimum_required(VERSION 3.5) +project(butter_worth_filter_design) + +set(CMAKE_TOOLCHAIN_FILE ${CMAKE_SOURCE_DIR}/toolchain-arm64.cmake) + + +include_directories(${CMAKE_SOURCE_DIR}) + +# 生产动态库 +add_library(butterworth STATIC + Butterworth.cpp + Biquad.cpp + conv.cpp + sos2tf.cpp + ) + +# 添加编译选项，按位置无关代码编译 +target_compile_options(butterworth PRIVATE -fPIC) + +# 添加测试 +enable_testing() +add_subdirectory(tests) diff --git a/README.md b/README.md index 3c411d4..feea772 100644 --- a/README.md +++ b/README.md @@ -10,7 +10,8 @@ The generated filter coefficients are split out into cascaded biquad sections, f * High order Parametric boost/cut EQ filter design *Biquad and Biquad Chain implementations (for filtering audio buffers with cascaded biquad sections) * Compact, readable and well-commented codebase - + * Convert from SOS to Transfer Function coefficients + * Convoluttion function ### Unit tests As with any good audio signal processing toolkit, there are unit tests that provide basic proof of correctness. There are currently 6 primary test cases that check 113 different conditions. @@ -19,9 +20,7 @@ Unit tests live in `main.cpp` and are written using the compact [Catch](https:// ### Prerequisites - * [SCONS](http://scons.org) as a cross-platform build system to build, test and run examples. - * On MacOS use [Homebrew](https://brew.sh): `$brew install scons` or [MacPorts](https://www.macports.org) `port install scons` - * On Linux: `apt-get install scons` + * cmake * [libsndfile](http://www.mega-nerd.com/libsndfile): `brew install libsndfile` ### Usage @@ -51,7 +50,24 @@ To generate the same set of coefficients in MATLAB (R14) as a comparison, to dou [sos, g] = zpk2sos(Zd, Pd, Kd) % zero-pole-gain form to second-order sections (SOS) ``` - +#### To convert from SOS to Transfer Function Form + +The Butterworth class get the Second Order Sections form coefficients of the Butter-Worth filter bank, but you can use `sos2tf()` to get the Transfer Function coefficients. + +``` +vector coeffs; // array of biquad filters (for this case, array size = 4 ) +Butterworth butterworth; +bool designedCorrectly = butterworth.loPass(44100, // fs + 500, // freq1 + 0, // freq2. N/A for lowpass + 8, // filter order, + coeffs, // coefficient array being filled + 1.0); // overall gain +// get Transfer Function coefficients +vectord b, a; +sos2tf(coeffs, gain, b, a); +// you can use b and a now. +``` ### Other filter design repos on GitHub * Vinnie Falco [DSP Filters](https://github.com/vinniefalco/DSPFilters) diff --git a/catch.hpp b/catch.hpp old mode 100755 new mode 100644 diff --git a/conv.cpp b/conv.cpp new file mode 100644 index 0000000..c4bb0fc --- /dev/null +++ b/conv.cpp @@ -0,0 +1,32 @@ + +#include "conv.hpp" + +/** + * Convolution of u and v vector. + */ +void conv(const vectord &u, const vectord &v, vectord &c) +{ + // fill with zeros + c.resize(u.size() + v.size() - 1, 0); + // use convolution machine implements convolution + for (int n = 0; n < c.size(); n++) + { + // iterate input signal u, kernal is v + for (int k = 0; k < v.size(); k++) + { + // just as the Math Equation of Convolution + int iu = n - k; + double uu; + if (iu < 0 || iu >= u.size()) + { + uu = 0; + } + else + { + uu = u[iu]; + } + // iterate kernel and cumulate it + c[n] = c[n] + v[k] * uu; + } + } +} diff --git a/conv.hpp b/conv.hpp new file mode 100644 index 0000000..0885522 --- /dev/null +++ b/conv.hpp @@ -0,0 +1,15 @@ +#ifndef CONV_H +#define CONV_H + +#include + +typedef std::vector vectord; + +/** + * Convolution of u and v vector. + * Return result vector c. + * length(c) = length(u) + length(v) - 1 + */ +void conv(const vectord& u, const vectord& v, vectord& c); + +#endif \ No newline at end of file diff --git a/main.cpp b/main.cpp index c3fdc17..15c909e 100644 --- a/main.cpp +++ b/main.cpp @@ -34,6 +34,7 @@ #define CATCH_CONFIG_MAIN // Tells Catch to provide a main() #include "catch.hpp" // Catch Unit test framework +#include using namespace std; diff --git a/sos2tf.cpp b/sos2tf.cpp new file mode 100644 index 0000000..88198cd --- /dev/null +++ b/sos2tf.cpp @@ -0,0 +1,67 @@ +#include "sos2tf.hpp" +#include "conv.hpp" + +// implement matlab statement: b(1:3) = sos(1, 1:3); +// exclude end index. +void slice_assign(vector &lhs, int lstart, int lend, vector &rhs, int rstart, int rend) +{ + for (int i = 0; i < lend - lstart; i++) + { + lhs[lstart + i] = rhs[rstart + i]; + } +} + +void sos2tf(const vector &sos, double gain, vector &b, vector &a) +{ + // Step1: Decompose b and a coefficients + int L = sos.size(); + + b.resize(2 * L + 1, 0); + a.resize(2 * L + 1, 0); + + // initializes the first polynomial coefficient + // b(1 : 3) = sos(1, 1 : 3); + // a(1 : 3) = sos(1, 4 : 6); + b[0] = sos[0].b0; + b[1] = sos[0].b1; + b[2] = sos[0].b2; + a[0] = 1; + a[1] = sos[0].a1; + a[2] = sos[0].a2; + + // iterate every Biquad + for (int r = 1; r < L; r++) + { + // get b and a of SOS + vectord b1(3), a1(3); + b1[0] = sos[r].b0; + b1[1] = sos[r].b1; + b1[2] = sos[r].b2; + a1[0] = 1; + a1[1] = sos[r].a1; + a1[2] = sos[r].a2; + // polynomial multiply by convolve. + // b(1:2*(r+1)+1) = myconv(b(1:2*r+1), b1); + std::vector::const_iterator b_begin = b.begin(); + std::vector::const_iterator b_end = b.begin() + 2 * r + 1; + vectord b0(b_begin, b_end); + vectord bc; + conv(b0, b1, bc); + slice_assign(b, 0, bc.size(), bc, 0, bc.size()); + + // a(1:2*(r+1)+1) = myconv(a(1:2*r+1), a1); + std::vector::const_iterator a_begin = a.begin(); + std::vector::const_iterator a_end = a.begin() + 2 * r + 1; + vectord a0(a_begin, a_end); + vectord ac; + conv(a0, a1, ac); + slice_assign(a, 0, ac.size(), ac, 0, ac.size()); + } + // cut size to 2L+1, because the convolution adds two zeros to the end + b.resize(2*L+1); + a.resize(2*L+1); + // multiply gain + for (int i = 0; i < b.size(); i++) { + b[i] *= gain; + } +} \ No newline at end of file diff --git a/sos2tf.hpp b/sos2tf.hpp new file mode 100644 index 0000000..5972914 --- /dev/null +++ b/sos2tf.hpp @@ -0,0 +1,45 @@ +#ifndef SOS2TF_H +#define SOS2TF_H + +#include "Biquad.h" + +/** + * My 2nd-order sections to transfer function model conversion. + * See matlab `sos2tf`. + * + * [B,A] = MYSOS2TF(SOS,GAIN) returns the numerator and denominator + * coefficients B and A of the discrete-time linear system given + * by the gain G and the matrix SOS in second-order sections form. + * + * SOS is an L by 6 matrix which contains the coefficients of each + * second-order section in its rows: + * SOS = [ b01 b11 b21 1 a11 a21 + * b02 b12 b22 1 a12 a22 + * ... + * b0L b1L b2L 1 a1L a2L ] + * + * The system transfer function is the product of the second-order + * transfer functions of the sections times the gain G. If G is not + * specified, it defaults to 1. Each row of the SOS matrix describes + * a 2nd order transfer function as + * b0k + b1k z^-1 + b2k z^-2 + * ---------------------------- + * 1 + a1k z^-1 + a2k z^-2 + * where k is the row index. + * + * In each row of the sos matrix, the first three elements form the + * coefficients of the numerator polynomial of a transfer function and + * the next three elements form the coefficients of the denominator + * polynomial. B and A are formed by multiplying all numerator and + * denominator polynomials respectively. Since each polynomial is of + * degree two, multiplying L such polynomials results in a polynomial of + * degree 2*L. + * + * Use convolution to multiply two polynomials. + * w = conv(u,v) Returns the convolution of the vectors u and v. + * If u and v are vectors of polynomial coefficients, convolving them is equivalent to multiplying the two polynomials. + * + */ +void sos2tf(const vector &sos, double gain, vector &b, vector &a); + +#endif \ No newline at end of file diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt new file mode 100644 index 0000000..ef5e573 --- /dev/null +++ b/tests/CMakeLists.txt @@ -0,0 +1,14 @@ +enable_testing() + +if(NOT ARM64) + add_executable(main-test main-test.cpp) + add_test(NAME main-test COMMAND main-test) + target_link_libraries(main-test PRIVATE butterworth -lsndfile) +endif() + +add_executable(test-butter test-butter.cpp) +add_test(NAME test-butter COMMAND test-butter) + +add_executable(test-conv test-conv.cpp) +add_test(NAME test-conv COMMAND test-conv) +target_link_libraries(test-conv PRIVATE butterworth) diff --git a/tests/main-test.cpp b/tests/main-test.cpp new file mode 100644 index 0000000..122c968 --- /dev/null +++ b/tests/main-test.cpp @@ -0,0 +1,873 @@ +/* + + This file is part of Butterworth Filter Design, a pair C++ classes and an + accompanying suite of unit tests for designing high order Butterworth IIR & + EQ filters using the bilinear transform. + The generated filter coefficients are split out into cascaded biquad sections, + for easy use in your garden variety biquad or second-order section (SOS). + + Reference: http://en.wikipedia.org/wiki/Butterworth_filter + http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt + + + Copyright (C) 2013, iroro orife + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . + + */ + + +#include +#include + +#include "Butterworth.h" +#include "sos2tf.hpp" +#include "conv.hpp" + +#define CATCH_CONFIG_MAIN // Tells Catch to provide a main() +#include "catch.hpp" // Catch Unit test framework +#include + +using namespace std; + + +TEST_CASE("Analogue lowpass prototype generation", " Verify poles match MATLAB buttap(8)"){ + + //****************************************************************************** + // + // baseline_poles values generated with MATLAB command >> buttap(8); + // + //****************************************************************************** + + int filterOrder = 8; + + vector poles; + poles.resize(filterOrder); + vector baseline_poles; + baseline_poles.resize(filterOrder); + + // Butterworth butter; + vector tempPoleStorage = Butterworth::prototypeAnalogLowPass(filterOrder); + + // copy tmppole into poles + int index = 0; + for(vector ::iterator itr = tempPoleStorage.begin(); itr != tempPoleStorage.end(); itr++){ + poles[index] = *itr; + index++; + } + + baseline_poles[0] = complex_double(-0.19509032201613, 0.98078528040323); + baseline_poles[1] = complex_double(-0.19509032201613, -0.98078528040323); + baseline_poles[2] = complex_double(-0.55557023301960, 0.83146961230255); + baseline_poles[3] = complex_double(-0.55557023301960, -0.83146961230255); + baseline_poles[4] = complex_double(-0.83146961230255, 0.55557023301960); + baseline_poles[5] = complex_double(-0.83146961230255, -0.55557023301960); + baseline_poles[6] = complex_double(-0.98078528040323, 0.19509032201613); + baseline_poles[7] = complex_double(-0.98078528040323, -0.19509032201613); + + // Tolerance + const double EPSILON = 1.0e-14; + + // Test + for(uint32_t i = 0; i < filterOrder; i++){ + CAPTURE(i); + REQUIRE(abs(baseline_poles[i].real() - poles[i].real()) <= EPSILON); + } +} + + +TEST_CASE("Lowpass filter coefficients", " Verify coefficients for a lowpass filter"){ + + //****************************************************************************** + // MATLAB (vR14) coefficients generated with the following code: + // + // [z, p, k] = butter(8, 500, 's'); + // [Zd, Pd, Kd] = bilinear(z, p, k, 44100); + // [sos, g] = zpk2sos(Zd, Pd, Kd) + //****************************************************************************** + + int filterOrder = 8; + double overallGain = 1.0; + + vector coeffs; // second-order sections (sos) + Butterworth butterworth; + + bool designedCorrectly = butterworth.loPass(44100, // fs + 500, // freq1 + 0, // freq2. N/A for lowpass + filterOrder, + coeffs, + overallGain); + REQUIRE(designedCorrectly == true); + + //****************************************************************************** + // MATLAB coefficients: first section + //****************************************************************************** + + double b0 = 1.0; + double b1 = 2.0; + double b2 = 1.0; + // double a0 = 1.0; + double a1 = -1.96762058043629 * (-1); // to convert to DF2T (direct form II transpose) + double a2 = 0.97261960500367 * (-1); // to convert to DF2T (direct form II transpose) + + int i = 0; + + const double EPSILON = 1.0e-4; + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); + + //****************************************************************************** + // MATLAB coefficients: second section + //****************************************************************************** + + b0 = 1.0; + b1 = 2.0; + b2 = 1.0; + // a0 = 1.0; + a1 = -1.91907529383978 * (-1); + a2 = 0.92395098208778 * (-1); + + i = 1; + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); + + //****************************************************************************** + // MATLAB coefficients: third section + //****************************************************************************** + + b0 = 1.0; + b1 = 2.0; + b2 = 1.0; + // a0 = 1.0; + a1 = -1.88350864178159 * (-1); + a2 = 0.88829396780773 * (-1); + + i = 2; + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); + + //****************************************************************************** + // MATLAB coefficients: fourth section + //****************************************************************************** + + b0 = 1.0; + b1 = 2.0; + b2 = 1.0; + // a0 = 1.0; + a1 = -1.86480445083537 * (-1); + a2 = 0.86954225616013 * (-1); + + i = 3; + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); +} + + +TEST_CASE("Highpass filter coefficients", " Verify coefficients for a highpass filter"){ + + //****************************************************************************** + // MATLAB (vR14) coefficients generated with the following code: + // + // [z, p, k] = butter(8, 500, 'high', 's'); + // [Zd, Pd, Kd] = bilinear(z, p, k, 44100); + // [sos, g] = zpk2sos(Zd, Pd, Kd) + //****************************************************************************** + + int filterOrder = 8; + double overallGain = 1.0; + const double EPSILON = 1.0e-4; + + vector coeffs; // second-order sections (sos) + Butterworth butterworth; + + bool designedCorrectly = butterworth.hiPass(44100, // fs + 500, // freq1 + 0, // freq2 + filterOrder, + coeffs, + overallGain); + REQUIRE(designedCorrectly == true); + + //****************************************************************************** + // MATLAB coefficients: first section + //****************************************************************************** + + double b0 = 1.0; + double b1 = -2.0; + double b2 = 1.0; + // double a0 = 1.0; + double a1 = -1.96762058043629 * (-1); + double a2 = 0.97261960500367 * (-1); + + int i = 0; + + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); + + //****************************************************************************** + // MATLAB coefficients: second section + //****************************************************************************** + + b0 = 1.0; + b1 = -2.0; + b2 = 1.0; + // a0 = 1.0; + a1 = -1.91907529383978 * (-1); + a2 = 0.92395098208778 * (-1); + + i = 1; + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); + + //****************************************************************************** + // MATLAB coefficients: third section + //****************************************************************************** + + b0 = 1.0; + b1 = -2.0; + b2 = 1.0; + // a0 = 1.0; + a1 = -1.88350864178159 * (-1); + a2 = 0.88829396780773 * (-1); + + i = 2; + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); + + //****************************************************************************** + // MATLAB coefficients: fourth section + //****************************************************************************** + + b0 = 1.0; + b1 = -2.0; + b2 = 1.0; + // a0 = 1.0; + a1 = -1.86480445083537 * (-1); + a2 = 0.86954225616013 * (-1); + + i = 3; + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); +} + + +TEST_CASE("Parametric EQ coefficients", " Verify coefficients for a parametric filter"){ + + //****************************************************************************** + // MATLAB coefficients were generated with the following code in + // MATLAB R14 with code hpeq1.m from + // + // Reference: Sophocles J. Orfanidis, "High-Order Digital Parametric Equalizer + // Design," J. Audio Eng. Soc., vol.53, pp. 1026-1046, November 2005. + //****************************************************************************** + + int filterOrder = 8; + double overallGain = 12.0; + const double EPSILON = 1.0e-4; + + vector coeffs; + Butterworth butterworth; + + bool designedCorrectly = butterworth.coefficientsEQ(Butterworth::kParametric, + 40000, // fs + 3500, // freq1 + 6500, // freq2 Df = 3kHz + filterOrder, + coeffs, + overallGain); + + REQUIRE(designedCorrectly == true); + + //****************************************************************************** + // MATLAB coefficients: first section + //****************************************************************************** + + double b0 = 1.0356; + double b1 = -2.6699; + double b2 = 3.4403; + double b3 = -2.3905; + double b4 = 0.8435; + // double a0 = 1.0; + double a1 = -2.6478 * (-1); + double a2 = 3.4810 * (-1); + double a3 = -2.4127 * (-1); + double a4 = 0.8384 * (-1); + + int i = 0; + + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(b3 - coeffs[i].b3) <= EPSILON); + REQUIRE(abs(b4 - coeffs[i].b4) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); + REQUIRE(abs(a3 - coeffs[i].a3) <= EPSILON); + REQUIRE(abs(a4 - coeffs[i].a4) <= EPSILON); + + //****************************************************************************** + // MATLAB coefficients: second section + //****************************************************************************** + + b0 = 1.0555; + b1 = -2.5476; + b2 = 2.9933; + b3 = -1.8553; + b4 = 0.5794; + // a0 = 1.0; + a1 = -2.4927 * (-1); + a2 = 3.0287 * (-1); + a3 = -1.9102 * (-1); + a4 = 0.5995 * (-1); + + i = 1; + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(b3 - coeffs[i].b3) <= EPSILON); + REQUIRE(abs(b4 - coeffs[i].b4) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); + REQUIRE(abs(a3 - coeffs[i].a3) <= EPSILON); + REQUIRE(abs(a4 - coeffs[i].a4) <= EPSILON); +} + + +TEST_CASE("Low shelf EQ coefficients", " Verify coefficients for a low shelf filter"){ + + //****************************************************************************** + // MATLAB coefficients were generated with the following code in + // MATLAB R14 with code hpeq1.m from + // + // Reference: Sophocles J. Orfanidis, "High-Order Digital Parametric Equalizer + // Design," J. Audio Eng. Soc., vol.53, pp. 1026-1046, November 2005. + //****************************************************************************** + + int filterOrder = 4; + double overallGain = 12.0; + const double EPSILON = 1.0e-4; + + vector coeffs; // second-order sections (sos) + Butterworth butterworth; + bool designedCorrectly = butterworth.coefficientsEQ(Butterworth::kLoShelf, + 40000, // fs + 0, // freq1 + 13000, // freq2 Df = 13kHz + filterOrder, + coeffs, + overallGain); + + REQUIRE(designedCorrectly == true); + + //****************************************************************************** + // MATLAB coefficients: first section + //****************************************************************************** + + double b0 = 1.6392; + double b1 = 1.7241; + double b2 = 0.9170; + // double a0 = 1.0; + double a1 = 0.6565 * (-1); + double a2 = 0.4887 * (-1); + + int i = 0; + + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); + + //****************************************************************************** + // MATLAB coefficients: second section + //****************************************************************************** + + b0 = 1.5790; + b1 = 1.2663; + b2 = 0.2984; + // a0 = 1.0; + a1 = 0.4822 * (-1); + a2 = 0.0934 * (-1); + + i = 1; + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); + + + //****************************************************************************** + // test another set of lowshelf design requirements + //****************************************************************************** + + designedCorrectly = butterworth.coefficientsEQ(Butterworth::kLoShelf, + 40000, // fs + 2000, // freq1 + 2800, // freq2 Df = 0.8kHz + filterOrder, + coeffs, + overallGain); + + REQUIRE(designedCorrectly == true); + + //****************************************************************************** + // MATLAB coefficients: first section + //****************************************************************************** + + b0 = 1.0222; + b1 = -1.8880; + b2 = 0.8949; + // a0 = 1.0; + a1 = -1.8953 * (-1); + a2 = 0.9099 * (-1); + + i = 0; + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); + + + //****************************************************************************** + // MATLAB coefficients: second section + //****************************************************************************** + + b0 = 1.0456; + b1 = -1.7749; + b2 = 0.7566; + // a0 = 1.0; + a1 = -1.7817 * (-1); + a2 = 0.7954 * (-1); + + i = 1; + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); +} + + +//------------------------------------------------------------------------------------------------- +TEST_CASE("High shelf EQ coefficients", " Verify coefficients for a high shelf filter"){ + + //****************************************************************************** + // MATLAB coefficients were generated with the following code in + // MATLAB R14 with code hpeq1.m from + // + // Reference: Sophocles J. Orfanidis, "High-Order Digital Parametric Equalizer + // Design," J. Audio Eng. Soc., vol.53, pp. 1026-1046, November 2005. + //****************************************************************************** + + int filterOrder = 4; + double overallGain = 12.0; + const double EPSILON = 1.0e-4; + + + vector coeffs; // second-order sections (sos) + Butterworth butterworth; + bool designedCorrectly = butterworth.coefficientsEQ(Butterworth::kHiShelf, + 40000, // fs + 0, // freq1 + 13000, // freq2 Df = 13kHz + filterOrder, + coeffs, + overallGain); + + + REQUIRE(designedCorrectly == true); + + //****************************************************************************** + // MATLAB coefficients: first section + //****************************************************************************** + + double b0 = 1.6392; + double b1 = -1.7241; + double b2 = 0.9170; + // double a0 = 1.0; + double a1 = -0.6565 * (-1); + double a2 = 0.4887 * (-1); + + int i = 0; + + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); + + //****************************************************************************** + // MATLAB coefficients: second section + //****************************************************************************** + + b0 = 1.5790; + b1 = -1.2663; + b2 = 0.2984; + // a0 = 1.0; + a1 = -0.4822 * (-1); + a2 = 0.0934 * (-1); + + i = 1; + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); + + + + //****************************************************************************** + // test another set of highshelf design requirements + //****************************************************************************** + + designedCorrectly = butterworth.coefficientsEQ(Butterworth::kHiShelf, + 40000, // fs + 2000, // freq1 + 2800, // freq2 Df = 13kHz + filterOrder, + coeffs, + overallGain); + + REQUIRE(designedCorrectly == true); + + //****************************************************************************** + // MATLAB coefficients: first section + //****************************************************************************** + + b0 = 1.0222; + b1 = 1.8880; + b2 = 0.8949; + // a0 = 1.0; + a1 = 1.8953 * (-1); + a2 = 0.9099 * (-1); + + i = 0; + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); + + //****************************************************************************** + // MATLAB coefficients: second section + //****************************************************************************** + + b0 = 1.0456; + b1 = 1.7749; + b2 = 0.7566; + // a0 = 1.0; + a1 = 1.7817 * (-1); + a2 = 0.7954 * (-1); + + i = 1; + REQUIRE(abs(b0 - coeffs[i].b0) <= EPSILON); + REQUIRE(abs(b1 - coeffs[i].b1) <= EPSILON); + REQUIRE(abs(b2 - coeffs[i].b2) <= EPSILON); + REQUIRE(abs(a1 - coeffs[i].a1) <= EPSILON); + REQUIRE(abs(a2 - coeffs[i].a2) <= EPSILON); +} + + +//------------------------------------------------------------------------------------------------- +#define BUFFER_LEN 1024 +#define MAX_CHANNELS 6 + +static void process_data (double *data, int count, int channels) +{ + int k; + + /* Process the data here. + ** If the soundfile contains more then 1 channel you need to take care of + ** the data interleaving youself. + ** Current we just apply a channel dependant gain. + */ + + for (k = 0 ; k < count ; k+= 1) + data [k] *= 1.0; + + return ; +} + +static void read_write_file (const char *infilename, const char *outfilename) +{ + static double data [BUFFER_LEN] ; + SNDFILE *infile, *outfile ; + SF_INFO sfinfo ; + int readcount ; + memset (&sfinfo, 0, sizeof (sfinfo)) ; + + if (! (infile = sf_open (infilename, SFM_READ, &sfinfo))) + { + printf ("Not able to open input file %s.\n", infilename) ; + puts (sf_strerror (NULL)) ; + return; + } + + if (sfinfo.channels > MAX_CHANNELS) + { printf ("Not able to process more than %d channels\n", MAX_CHANNELS) ; + return; + } + if (! (outfile = sf_open (outfilename, SFM_WRITE, &sfinfo))) + { printf ("Not able to open output file %s.\n", outfilename) ; + puts (sf_strerror (NULL)) ; + return; + } + + while ((readcount = sf_read_double (infile, data, BUFFER_LEN))) + { process_data (data, readcount, sfinfo.channels) ; + sf_write_double (outfile, data, readcount) ; + } + + sf_close (infile); sf_close (outfile) ; + return; +} + + +TEST_CASE("Test gain control issue #3"){ + + //****************************************************************************** + // For this test you will need libsndfile.a which lives in + // /opt/local/lib/libsndfile.a installed on MacOS via macports + // https://trac.macports.org/browser/trunk/dports/audio/libsndfile/Portfile + //****************************************************************************** + + const char * infilename = "../../sweep_0_20K.wav" ; + const char * outfilename = "sweep_0_20K_filtered.wav" ; + read_write_file(infilename,outfilename); +} + +TEST_CASE("Bandpass test") { + + int filterOrder = 4; + const double EPSILON = 1.0e-4; + double gain; + + vector coeffs; // second-order sections (sos) + Butterworth butterworth; + bool designedCorrectly = butterworth.bandPass( + 100, // fs + 10, // freq1 + 0.5, // freq2 + filterOrder, + coeffs, gain); + + cout << "gain: " << gain << endl; + for(int i=0; i<4; i++) { + cout << " === section " << i << " ===" << endl; + cout << "\tBiquad parameters b0: " << coeffs[i].b0 << endl; + cout << "\tBiquad parameters b1: " << coeffs[i].b1 << endl; + cout << "\tBiquad parameters b2: " << coeffs[i].b2 << endl; + cout << "\tBiquad parameters a1: " << coeffs[i].a1 << endl; + cout << "\tBiquad parameters a2: " << coeffs[i].a2 << endl; + } +} + +TEST_CASE("Convolution u and v", " Verify convolution is correct") +{ + const double EPSILON = 1E-4; + + vectord u = {1, 2, 3, 4, 5}; + vectord v = {1, 2, 3}; + vectord c; + conv(u, v, c); + + vectord r = {1, 4, 10, 16, 22, 22, 15}; + for (size_t i = 0; i < c.size(); i++) + { + REQUIRE(abs(r[i] - c[i]) <= EPSILON); + } +} + +TEST_CASE("SOS To Transfer Function coefficients", " Verify coefficients is correct") +{ + const double EPSILON = 1E-4; + + vector sos(4); + sos[0] = Biquad(1, 2, 1, -1.11879834477625, 0.336688180692519); + sos[1] = Biquad(1, 2, 1, -1.35028339447298, 0.660543824885074); + sos[2] = Biquad(1, -2, 1, -1.93859735792819, 0.939736057085694); + sos[3] = Biquad(1, -2, 1, -1.97691413926085, 0.977914024909786); + + double gain = 0.00404923315402386; + + vectord b, a; + sos2tf(sos, gain, b, a); + + // matlab b,a + vectord mb = {0.00404923315402386, 0, -0.0161969326160954, 0, 0.0242953989241431, 0, -0.0161969326160954, 0, 0.00404923315402386}; + vectord ma = {1, -6.38459323643827, 17.9257353790401, -28.9644589090336, 29.5037393371934, -19.4170269899449, 8.06393935848502, -1.93171376981280, 0.204378907481348}; + + for (int i = 0; i < b.size(); i++) + { + CAPTURE(i) + REQUIRE(abs(b[i] - mb[i]) <= EPSILON); + REQUIRE(abs(a[i] - ma[i]) <= EPSILON); + } +} + +TEST_CASE("100Hz Butter-Worth Band-Pass Filter Design Output Transfer Function Coefficients", "Verify coefficients with matlab values") +{ + int filterOrder = 4; + double gain = 1.0; + double Fs = 100; + double fp1 = 0.5; + double fs1 = 10; + + vector coeffs; // second-order sections (sos) + Butterworth butterworth; + + bool designedCorrectly = butterworth.bandPass(Fs, // fs + fp1, // freq1 + fs1, // freq2. N/A for lowpass + filterOrder, + coeffs, + gain); + REQUIRE(designedCorrectly == true); + + // Inverse sign of a1, a2 + // Because we dont know why Butterworth class always gets negative of a1 and a2 + for (int i=0; i coeffs; // second-order sections (sos) + Butterworth butterworth; + + bool designedCorrectly = butterworth.bandPass(Fs, // fs + fp1, // freq1 + fs1, // freq2. N/A for lowpass + filterOrder, + coeffs, + gain); + REQUIRE(designedCorrectly == true); + + // Inverse sign of a1, a2 + // Because we dont know why Butterworth class always gets negative of a1 and a2 + for (int i=0; i coeffs; // second-order sections (sos) + Butterworth butterworth; + + bool designedCorrectly = butterworth.loPass(Fs, // fs + fp1, // freq1 + fs1, // freq2. N/A for lowpass + filterOrder, + coeffs, + gain); + REQUIRE(designedCorrectly == true); + + // Inverse sign of a1, a2 + // Because we dont know why Butterworth class always gets negative of a1 and a2 + for (int i=0; i + + +using namespace std; + +int main() +{ + cout << "Butter-Worth Filter: " << endl; +} \ No newline at end of file diff --git a/tests/test-conv.cpp b/tests/test-conv.cpp new file mode 100644 index 0000000..55e465e --- /dev/null +++ b/tests/test-conv.cpp @@ -0,0 +1,48 @@ +#define CATCH_CONFIG_MAIN // Tells Catch to provide a main() +#include "catch.hpp" // Catch Unit test framework +#include +#include "conv.hpp" +#include "sos2tf.hpp" + +using namespace std; + +const double EPSILON = 1E-4; + +TEST_CASE("Convolution u and v", " Verify convolution is correct") +{ + vectord u = {1, 2, 3, 4, 5}; + vectord v = {1, 2, 3}; + vectord c; + conv(u, v, c); + + vectord r = {1, 4, 10, 16, 22, 22, 15}; + for (size_t i = 0; i < c.size(); i++) + { + REQUIRE(abs(r[i] - c[i]) <= EPSILON); + } +} + +TEST_CASE("SOS To Transfer Function coefficients", " Verify coefficients is correct") +{ + vector sos(4); + sos[0] = Biquad(1, 2, 1, -1.11879834477625, 0.336688180692519); + sos[1] = Biquad(1, 2, 1, -1.35028339447298, 0.660543824885074); + sos[2] = Biquad(1, -2, 1, -1.93859735792819, 0.939736057085694); + sos[3] = Biquad(1, -2, 1, -1.97691413926085, 0.977914024909786); + + double gain = 0.00404923315402386; + + vectord b, a; + sos2tf(sos, gain, b, a); + + // matlab b,a + vectord mb = {0.00404923315402386, 0, -0.0161969326160954, 0, 0.0242953989241431, 0, -0.0161969326160954, 0, 0.00404923315402386}; + vectord ma = {1, -6.38459323643827, 17.9257353790401, -28.9644589090336, 29.5037393371934, -19.4170269899449, 8.06393935848502, -1.93171376981280, 0.204378907481348}; + + for (int i = 0; i < b.size(); i++) + { + CAPTURE(i) + REQUIRE(abs(b[i] - mb[i]) <= EPSILON); + REQUIRE(abs(a[i] - ma[i]) <= EPSILON); + } +} diff --git a/toolchain-arm64.cmake b/toolchain-arm64.cmake new file mode 100644 index 0000000..92004c7 --- /dev/null +++ b/toolchain-arm64.cmake @@ -0,0 +1,23 @@ +set(CMAKE_SYSTEM_NAME Linux) # Operating system type +set(CMAKE_SYSTEM_PROCESSOR aarch64) # CPU architecture +set(CMAKE_C_COMPILER_TARGET aarch64-linux-gnu) + +set(CMAKE_C_COMPILER aarch64-linux-gnu-gcc) # gcc compiler program name +set(CMAKE_CXX_COMPILER aarch64-linux-gnu-g++) # g++ compiler program name +set(CMAKE_Fortran_COMPILER aarch64-linux-gnu-gfortran) # g++ compiler program name + +set(CMAKE_FIND_ROOT_PATH /opt/cross-gcc-linaro-6.3.1-arm64/) # Root directory of g++ +set(CMAKE_FIND_ROOT_PATH_MODE_PROGRAM BOTH) # Tell CMake to look for programs like g++ from under ROOT_PATH, but it can also look elsewhere +set(CMAKE_FIND_ROOT_PATH_MODE_LIBRARY ONLY) # Tell CMake to look for library files from under ROOT_PATH +set(CMAKE_FIND_ROOT_PATH_MODE_INCLUDE ONLY) # Tell CMake to look for header files from under ROOT_PATH +# set(CMAKE_PREFIX_PATH /opt/cross-gcc-linaro-6.3.1-arm64/arm-linux-gnueabihf) + +# Whether to statically link libc and libstdc++ +# set(CMAKE_C_FLAGS "-static-libgcc -static-libstdc++") + +# Set compilation options +# set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -march=armv7-a -mfpu=vfpv3") +# set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=armv7-a -mfpu=vfpv3") # fpu版本 + +# Set link options +set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -Wl,-rpath-link,/opt/cross-gcc-linaro-6.3.1-arm64/aarch64-linux-gnu/lib") diff --git a/toolchain-armv7-32b.cmake b/toolchain-armv7-32b.cmake new file mode 100644 index 0000000..f18fd11 --- /dev/null +++ b/toolchain-armv7-32b.cmake @@ -0,0 +1,22 @@ +set(CMAKE_SYSTEM_NAME Linux) # Operating system type +set(CMAKE_SYSTEM_PROCESSOR arm) # CPU architecture + +set(CMAKE_C_COMPILER arm-linux-musleabihf-gcc) # gcc compiler program name +set(CMAKE_CXX_COMPILER arm-linux-musleabihf-g++) # g++ compiler program name + +set(CMAKE_FIND_ROOT_PATH /opt/cross/) # Root directory of g++ +set(CMAKE_FIND_ROOT_PATH_MODE_PROGRAM BOTH) # Tell CMake to look for programs like g++ from under ROOT_PATH, but it can also look elsewhere +set(CMAKE_FIND_ROOT_PATH_MODE_LIBRARY ONLY) # Tell CMake to look for library files from under ROOT_PATH +set(CMAKE_FIND_ROOT_PATH_MODE_INCLUDE ONLY) # Tell CMake to look for header files from under ROOT_PATH + +set(CMAKE_SYSTEM_NAME Linux) +set(CMAKE_C_COMPILER_TARGET arm-linux-musleabihf) # Target programs that use the musl libc library, if they are standard gnu libraries, can comment out this line + +# set(CMAKE_C_FLAGS "-static-libgcc -static-libstdc++") + +# Set compilation options +set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -march=armv7-a -mfpu=vfpv3") +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=armv7-a -mfpu=vfpv3") # fpu版本 + +# Set link options +set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} -Wl,-rpath-link,/opt/cross/arm-linux-musleabihf/lib")