Matlab fdatool设计实现一维数字滤波器

主旨

Matlab的fdatool生成的滤波器系数与滤波器传输函数的A/B参数的对应关系并不直观，本文目的在于总结出滤波器系数到传输函数A/B参数的转换方法和原理。

问题描述

fdatool是Matlab提供的一个交互式数字滤波器设计工具[1]. 这个工具可以根据用户给定的滤波器类型、通带阻带参数等信息设计滤波器，并将滤波器的各类响应曲线(幅频响应、相频响应、冲击响应、延迟等)画出来。

当得到响应符合预期的滤波器设计之后，下一步需要将滤波器参数导出并代入实际运算。Matlab的一维数字滤波器计算函数为filter（无时延补偿）或filtfilt(有时延补偿)，其需要的参数除了原始信号外，还包括传输函数的A/B参数。其中A/B参数分别代表如下的滤波器传输函数表达式中数列A和数列B

H(z)=b(1)+b(2)z−1+..b(nb+1)z−nba(1)+a(2)z−1+..a(na+1)z−na

FDATool生成的系数并不能直观的对应于上述A/B参数。对于FIR滤波器，Export得到的是一个1*(N+1)的向量；对于IIR滤波器，export得到的是两个矩阵，一个L*6, 另一个为(L+1)*1. 如何将上述参数转换为A/B参数就是本文要讨论和解决的问题

FIR滤波器

FIR滤波器基本结构如下图所示[3]




By BlanchardJ - Own work, Public Domain, Link

上图中不存在从y(n)的反馈，所以其传输函数可以简化为

H(z)=b(1)+b(2)z−1+..b(nb+1)z−nb

因此，FdaTool生成的1*(N+1)向量可以对应于上述B参数，而A参数只要输入1*1向量，且a(1) = 1即可。Matlab参考代码

b = Num; %Num为fdatool输出的向量
a =  [1];
filteredData = filtfilt(b, a, x); %x为滤波器输入信号

IIR滤波器

IIR滤波器结构如下图所示[4]




CC BY-SA 3.0, Link

这里就存在反馈，因此参数A不再是单一元素。针对IIR滤波器，Fdatool输出的sos向量(L*6)是以second order section representation表示的滤波器系数，而g向量((L+1)*1)是增益向量。

所谓second order section，指的是以一系列级联的的biquadratric section表示的IIR滤波器[5]

A common strategy is to realize a higher-order (greater than 2) digital filter as a cascaded series of second-order “biquadratric” (or “biquad”) sections[6] (see digital biquad filter).

而biquadratric section本身是一种二阶递归线性滤波器[6]，其传输方程如下

H(z)=b(1)+b(2)z−1+b(3)z−21+a(2)z−1+a(2)z−2

整个IIR滤波器的传输函数可以写为如下形式，这里g(L+1)=1，因此没有写入下面的公式

H(z)=∏i=1Lg(i)b(i,1)+b(i,2)z−1+b(i,3)z−21+a(i,2)z−1+a(i,3)z−2

如何将上述形式的sos向量和g向量转化为滤波器传输函数的一般形式？Matlab提供了sos2tf功能。样例代码如下

[b, a] = sos2tf(SOS,G);%将sos向量和g向量转换为A/B参数
filteredData = filtfilt(b, a, x); %x为滤波器输入信号

参考资料

[1] https://cn.mathworks.com/help/signal/examples/introduction-to-filter-designer.html

[2] http://www.mathworks.com/help/matlab/ref/filter.html

[3] https://en.wikipedia.org/wiki/Finite_impulse_response

[4] https://en.wikipedia.org/wiki/Infinite_impulse_response

[5] https://en.wikipedia.org/wiki/Digital_filter#Cascaded_second-order_sections

[6] https://en.wikipedia.org/wiki/Digital_biquad_filter