声明：感谢知乎大佬的文章，原文链接

数字滤波器实现方法是把滤波器所要完成的运算编成程序并让计算机执行,也就是采用在代码的形式。
它面对的是离散时间的数字信号，是把输入序列通过一定的运算变换成输出序列。
问：如何将连续的模拟滤波器变成离散的数字滤波器？
答：双线性变换
$$S=\frac{2}{Ts}\frac{1-z^{-1}}{1+z^{-1}}=2fs\frac{1-z^{-1}}{1+z^{-1}}$$
S域：复频域 Z域：极坐标系

一阶数字滤波器

时域分析

电容
$$C=\frac{Q}{V}（单位电压下器件所存储的电荷量）$$
瞬时电流
$$Ic=\frac{dQ}{dt}（单位时间内通过的电荷数）$$
代入得
$$Ic=\frac{d(C\ast V)}{dt}$$
化简得
$$Ic=C\ast \frac{d{U}_o}{dt}$$
根据基尔霍夫定律
$$Ui=Ic\ast R+U_o$$
代入化简得
$$Ui=RC\ast \frac{d{U}_o}{dt}+U_o$$

$$\frac{d{U}_o}{dt}+\frac{U_o}{RC}=\frac{U_i}{RC}$$
根据微分公式
$$y^{\prime }+p(x)y=q(x)$$

$$y=e^{-\int p(x)dx}[\int e^{\int p(x)dx}q(x)+C]$$

$$U_o=e^{-\int \frac{1}{RC}dt}[\int e^{\int \frac{1}{RC}dt}\frac{U_i}{RC}+C]$$

$$U_o=e^{-\frac{t}{RC}}[\frac{U_i}{RC}\int e^{\frac{t}{RC}dt}+C]$$
其中
$$\int e^{t}dt=e^t$$

$$U_o=e^{-\frac{t}{RC}}[\frac{Ui}{RC}\int RC\ast e^{\frac{t}{RC}d\frac{t}{RC}}+C]$$

$$U_o=e^{-\frac{t}{RC}}[\frac{Ui}{RC}\ast RC\ast e^{\frac{t}{RC}}+C]$$

$$U_o=C\ast e^{-\frac{t}{RC}}+U_i$$

$$当t=0时，U_o=0,得出C=-U_i$$

$$U_o=（1-e^{-\frac{t}{RC}}）\ast U_i$$

频域分析

电容的阻抗
$$R=\frac{1}{jwC}$$
电压的传递函数
$$\frac{U_o}{U_i}=\frac{\frac{1}{jwC}}{R+\frac{1}{jwC}}$$
$$\frac{U_o}{U_i}=\frac{1}{1+jwRC}$$
化简
$$\frac{U_o}{U_i}=\frac{1}{1+(wRC{)}^2}-j\frac{wRC}{1+(wRC{)}^2}$$
对其取模得
$$\sqrt{(\frac{1}{1+(wRC{)}^2}{)}^2+(\frac{wRC}{1+(wRC{)}^2}{)}^2}$$
$$\sqrt{\frac{1}{1+(wRC{)}^2}}$$
$$当\frac{Uo}{Ui}=\frac{\sqrt{2}}{2}，此时的频率为截止频率$$
$$\frac{\sqrt{2}}{2}=\sqrt{\frac{1}{1+(wRC{)}^2}}$$
$$wRC=1,w=2\pi f,得出2\pi fRC=1$$
截止频率
$$f=\frac{1}{2\pi RC}$$

数字化

$$H(jw)=\frac{U_o}{U_i}=\frac{1}{1+jwRC}$$
$$s=jw$$
$$H(s)=\frac{1}{1+RCs}$$
离散化 一阶后向差分进入Z变换
$$s=\frac{1-z^{-1}}{T}（T:采样周期）$$
$$H(z)=\frac{U_o(z)}{U_i(z)}=\frac{1}{1+RC\frac{1-z^{-1}}{T}}$$
化简得
$$H(z)=\frac{U_o(z)}{U_i(z)}=\frac{1}{1+RC\frac{1-z^{-1}}{T}}$$
$$U_o(z)=\frac{T}{RC+T}{U}_i(z)+\frac{RC}{RC+T}{U}_o(z)z^{-1}$$
Z反变换得
$$U_o(n)=\frac{T}{RC+T}{U}_i(n)+\frac{RC}{RC+T}{U}_o(n-1)$$
变换得
$$U_o(n)=\frac{T}{RC+T}{U}_i(n)+(1-\frac{T}{RC+T})U_o(n-1)$$
$$令A=\frac{T}{RC+T}，U_o(n)=A{U}_i(n)+(1-A)U_o(n-1)$$
前面的频域得出
$$f=\frac{1}{2\pi RC},求得RC=\frac{1}{2\pi f}$$
从而得出
$$A=\frac{1}{1+\frac{1}{2\pi fT}}$$
$$此时将A代入Uo(n)中,得到输入Ui在截止频率f输出Uo在数字域中的一阶滤波器表达式$$
$$U_o(n)=\frac{1}{1+\frac{1}{2\pi fT}}{U}_i(n)+(1-\frac{1}{1+\frac{1}{2\pi fT}})U_o(n-1)$$
$$U_o(n)=\frac{1}{1+\frac{1}{2\pi fT}}(U_i(n)-U_o(n-1))+U_o(n-1)$$

代码示例

截止频率 f = 1Hz
thr_lpf+=(1 / (1 + 1/(2.0f * 3.14f * T )))*(height_thr - thr_lpf) ，其中T(采样周期)=1/采样频率

二阶巴特沃斯低通滤波器

S域和Z域的频率关系分析

对于s域来说，在模拟截止频率为 fa 时有：
$$S=j{w}_a=j2\pi f_a$$
对于z域来说，在数字截止频率为 fd 的情况下
$$z=e^{sT},其中s=j{w}_d,w_d=2\pi f_d,T={\frac{1}{f}}_s$$
$$z=e^{j2\pi \frac{f_d}{f_s}}$$
双线性变换
$$S=\frac{2}{T_s}\frac{1-z^{-1}}{1+z^{-1}}=2f_s\frac{1-z^{-1}}{1+z^{-1}}$$
把这个数字截止频率为 fd 带入双线性变换公式可以得到：
$$S=2f_s\frac{1-z^{-1}}{1+z^{-1}}=2f_s\frac{1-e^{-j2\pi \frac{f_d}{f_s}}}{1+e^{-j2\pi \frac{f_d}{f_s}}}=2f_s\frac{e^{j2\pi \frac{f_d}{f_s}}-1}{e^{j2\pi \frac{f_d}{f_s}}+1}$$
运用下方的公式进行化简
$$e^{jx}=cosx+jsinx$$
$$\frac{e^{jA}-e^{-jA}}{e^{jA}+e^{-jA}}=\frac{cosA+jsinA-(cos(-A)+jsin(-A))}{cosA+jsinA+cos(-A)+jsin(-A)}=jtanA$$
$$\frac{e^{jw}-1}{e^{jw}+1}=\frac{e^{\frac{jw}{2}}}{e^{\frac{jw}{2}}}\ast \frac{e^{\frac{jw}{2}}-e^{-\frac{jw}{2}}}{e^{\frac{jw}{2}}+e^{-\frac{jw}{2}}}=\frac{e^{\frac{jw}{2}}-e^{-\frac{jw}{2}}}{e^{\frac{jw}{2}}+e^{-\frac{jw}{2}}}=jtan\frac{w}{2}$$
化简得
$$S=2f_s\frac{e^{j2\pi \frac{f_d}{f_s}}-1}{e^{j2\pi \frac{f_d}{f_s}}+1}=j2f_{s}tan\frac{\pi f_d}{f_s}$$
$$S=j2\pi f_a=j2f_{s}tan\frac{\pi f_d}{f_s}$$
$$f_a=\frac{f_s}{\pi }tan\frac{\pi f_d}{f_s}$$
模拟滤波器的角频率Wa
$$w_a=2\pi f_a=2\pi \ast \frac{f_s}{\pi }tan\frac{\pi f_d}{f_s}=2f_{s}tan\frac{\pi f_d}{f_s}$$
去归一化只需要将模拟滤波器传递函数中的s进行如下替换
$$S=\frac{S}{w_a}$$
得出用数字截止频率 fd 表示S的公式
$$S=\frac{2f_s\frac{1-z^{-1}}{1+z^{-1}}}{2f_{s}tan\frac{\pi f_d}{f_s}}=\frac{\frac{1-z^{-1}}{1+z^{-1}}}{tan\frac{\pi f_d}{f_s}}$$

巴特沃斯滤波器举例说明

$$上表为巴特沃斯滤波器归一化（截止频率为1弧度\frac{1}{2}\pi )系数表，我们采用二阶(N=2)进行举例$$
$$H(s)=\frac{1}{S^2+1.414S+1}$$
$$H(z)=\frac{1}{(\frac{\frac{1-z^{-1}}{1+z^{-1}}}{tan\frac{\pi f_d}{f_s}}{)}^2+1.414(\frac{\frac{1-z^{-1}}{1+z^{-1}}}{tan\frac{\pi f_d}{f_s}})+1}$$
为方便化简将
$$Q=tan\frac{\pi f_d}{f_s}$$
化简得
$$H(z)=\frac{Y(z)}{X(z)}=\frac{Q^2+2Q^2{z}^{-1}+Q^2{z}^{-2}}{(1+1.414Q+Q^2)+(2Q^{-2}-2)z^{-1}+(1+Q^2-1.414Q)z^{-2}}$$
再次为方便化简将
$$B_0=Q^2,B_1=2Q^2,B_2=Q^2,C=1+1.414Q+Q^2,A_1=2Q^{-2}-2,A_2=1+Q^2-1.414Q$$
化简得
$$H(z)=\frac{Y(z)}{X(z)}=\frac{B_0+B_1{z}^{-1}+B_2{z}^{-2}}{C+A_1{z}^{-1}+A_2{z}^{-2}}$$
$$Y(z)C+A_{1}Y(z)z^{-1}+A_{2}Y(z)z^{-2}=B_{0}X(z)+B_{1}X(z)z^{-1}+B_{2}X(z)z^{-2}$$
$$Y(z)=\frac{B_0}{C}X(z)+\frac{B_1}{C}X(z)z^{-1}+\frac{B_2}{C}X(z)z^{-2}-\frac{A_1}{C}Y(z)z^{-1}-\frac{A_2}{C}Y(z)z^{-2}$$
再再次为方便化简将
$$b_0=\frac{B_0}{C},b_1=\frac{B_1}{C},b_2=\frac{B_2}{C},a_1=\frac{A_1}{C},a_2=\frac{A_2}{C}$$
$$Y(z)=b_{0}X(z)+b_{1}X(z^{-1})+b_{2}X(z^{-2})-a_{1}Y(z^{-1})-a2Y(z^{-2})$$
$$z:表示当前周期，z^{-1}:表示上个周期，z^{-2}:表示上上个周期$$

代码示例

//二阶滤波器系数结构体
struct DoubleFilterFactor
{float b0;float b1;float b2;float a1;float a2;
}
//二阶滤波器暂存X,Y数据结构体
struct XYData
{float XData[2];float YData[2];
}void double IIRDoubleFilter(struct DoubleFilterFactor* facotr, struct XYData* xydata, float in)
{float y = facotr->b0*in + facotr->b1*xydata->XData[0] +facotr->b2*xydata->XData[1] +facotr->a1*xydata->YData[0] +facotr->a2*xydata->YData[0];//y表示输出信号xydata->XData[1] = xydata->XData[0];xydata->XData[0] = in;xydata->YData[1] = xydata->YData[0];xydata->YData[0] = y;return y;
}使用方式：
struct DoubleFilterFactor* tmpfactor;//对该二阶滤波器系数结构体赋值
struct XYData* tmpxydata;//无需赋值
int i;
float sample[num];//x表示输入信号
for(i = 0;i < num;i++)
{double Filterdata = IIRDoubleFilter(tmpfactor, tmpxydata, sample[i]);//Filterdata就是经过二阶滤波后的数据
}