Added sos2tf() function to convert SOS to Transfer Function coefficients like matlab's SOS2TF

.gitignore

@@ -0,0 +1,3 @@
+/build/
+/.vscode/
+/build-*/

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(){
 }
 

Biquad.h

@@ -27,6 +27,8 @@
 
  */
 
+#ifndef BIQUAD_H
+#define BIQUAD_H
 
 #include <vector>
 #include <complex>
@@ -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

Butterworth.cpp

@@ -28,7 +28,7 @@
  */
 
 
-#import <iostream>
+#include <iostream>
 #include <math.h>
 #include "Butterworth.h"
 

Butterworth.h

@@ -27,6 +27,8 @@
 
  */
 
+#ifndef BUTTERWORTH_H
+#define BUTTERWORTH_H
 
 #include "Biquad.h"
 
@@ -116,3 +118,4 @@ class Butterworth {
     double * ba;
 };
 
+#endif

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)

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 <Biquad> 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)

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;
+        }
+    }
+}

conv.hpp

@@ -0,0 +1,15 @@
+#ifndef CONV_H
+#define CONV_H
+
+#include <vector>
+
+typedef std::vector<double> 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

main.cpp

@@ -34,6 +34,7 @@
 
 #define CATCH_CONFIG_MAIN   // Tells Catch to provide a main()
 #include "catch.hpp"        // Catch Unit test framework
+#include <cstring>
 
 using namespace std;
 

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<double> &lhs, int lstart, int lend, vector<double> &rhs, int rstart, int rend)
+{
+    for (int i = 0; i < lend - lstart; i++)
+    {
+        lhs[lstart + i] = rhs[rstart + i];
+    }
+}
+
+void sos2tf(const vector<Biquad> &sos, double gain, vector<double> &b, vector<double> &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<double>::const_iterator b_begin = b.begin();
+        std::vector<double>::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<double>::const_iterator a_begin = a.begin();
+        std::vector<double>::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;
+    }
+}

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<Biquad> &sos, double gain, vector<double> &b, vector<double> &a);
+
+#endif

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)