巴特沃斯一阶低通滤波器参数推导

1. 关于欧拉公式

1. 欧拉公式中的e就是自然对数的底数e
2. 欧拉公式不是仅仅是一种自定义的表示方法，是可以推导的
   首先将以下三个表达式用泰勒级数展开：
   $$e^x=1+x+\frac{x^2}{2!}+\frac{x^3}{3!}+...\text{\hspace{1em}(1)}$$

$$\cos (x)=1-\frac{x^2}{2!}+\frac{x^4}{4!}-\frac{x^6}{6!}...\text{\hspace{1em}(2)}$$

$$\sin (x)=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}...\text{\hspace{1em}(3)}$$

所以：
$$\begin{gathered} e^{ix}=1+ix+\frac{(ix{)}^2}{2!}+\frac{(ix{)}^3}{3!}+\frac{(ix{)}^4}{4!}+\frac{(ix{)}^5}{5!}+\frac{(ix{)}^6}{6!}+\frac{(ix{)}^7}{7!}... \\ =1+ix-\frac{x^2}{2!}-\frac{i{x}^3}{3!}+\frac{(x{)}^4}{4!}+\frac{i{x}^5}{5!}-\frac{x^6}{6!}-\frac{i{x}^7}{7!}... \\ =[1-\frac{x^2}{2!}+\frac{(x{)}^4}{4!}-\frac{x^6}{6!}...]+i[x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}...] \\ =\cos (x)+i\sin (x) \end{gathered}$$

2. 归一化频率

信号处理工具箱中经常使用的频率是Nyquist频率，它被定义为采样频率的一半，在滤波器的结束选择和设计当中的截止频率均使用Nyquist频率进行归一化处理。
例如，对于一个采样频率为1000Hz的系统，300Hz的归一化即为300/500=0.6。归一化频率的范围在[0,1]之间。如果要将归一化频率转换为角频率，则将归一化频率乘以pi；如果将归一化频率转换成Hz，则将归一化频率乘以采样频率的一半。
采样率的一半是最高频率，认为是１，那么真实频率和最高频率的比值就是归一化频率！也叫数字频率！

3. 巴特沃斯一阶低通滤波器参数推导

参考文献:飞控中的IIR二阶滤波器

z域与s域截止频率的关系

把一个传递函数从s域变换到z域时，如果变换方法是有差的（一般都是有差的，因为损失了时间分辨率），则变换前后的滤波特性（包括截止频率）会发生变化，此处以双线性变换为例。

双线性变换的表达式：
$$s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}}\text{\hspace{1em}(3.1)}$$

在上述变换中，再将z域表达式无损地变换回s域，（为了便于区分，给s或者$\omega$加上a和d下标分别表示初始系统和变换后的系统）：
$$z=e^{s_{d}T}=e^{j{\omega }_{d}T}\text{\hspace{1em}(3.2)}$$

将(3.2)代入(3.1)：
$$s_a=j{\omega }_a=2f_s\frac{1-e^{-j{\omega }_{d}T}}{1+e^{-j{\omega }_{d}T}}=2f_s\frac{e^{j{\omega }_{d}T}-1}{e^{j{\omega }_{d}T}+1}=j\frac{2}{T}\tan ({\omega }_{d}T/2)\text{\hspace{1em}(3.3)}$$

上式的推导过程：
根据欧拉公式：
$$e^{ix}=cosx+isinx\text{\hspace{1em}(3.4)}$$
可得：
$$\begin{gathered} \frac{e^{iA}-e^{-jA}}{e^{iA}+e^{-iA}}=\frac{cos(A)+isinA-(cos(-A)+isin(-A))}{cosA+isinA+cos(-A)+isin(-A)} \\ =\frac{2isinA}{2cosA}=itanA\text{\hspace{1em}(3.5)} \end{gathered}$$
而
$$\frac{e^{j{\omega }_{d}T}-1}{e^{j{\omega }_{d}T}+1}=\frac{e^{j{\omega }_{d}T/2}}{e^{j{\omega }_{d}T/2}}\frac{e^{j{\omega }_{d}T/2}-e^{-j{\omega }_{d}T/2}}{e^{j{\omega }_{d}T/2}+e^{-j{\omega }_{d}T/2}}=jtan({\omega }_{d}T/2)\text{\hspace{1em}(3.6)}$$

对(3.3)式两边取模可得：
$${\omega }_a=\frac{2}{T}\tan ({\omega }_{d}T/2)\text{\hspace{1em}(3.7)}$$

$$2\pi f_a=2f_s\tan (\frac{2\pi f_d}{2f_s})\Rightarrow f_a=\frac{f_s}{\pi }\tan (\frac{\pi f_d}{f_s})\text{\hspace{1em}(3.8)}$$
$f_d(x)$和$f_a(y)$的函数关系是这样的($f_s=4K$)：

可见频率越小，$f_a$和$f_d$越接近，频率越高，差别越大。

---

巴特沃斯低通滤波器的s域标准形式：
$$H_{an}(s)=\frac{d_0}{a_0+a_{1}s+a_2{s}^2+...+a_N{s}^N}\text{\hspace{1em}(3.9)}$$
一般取直流增益为1，此时有
$$d_0=a_0=1$$

对应的归一化参数表：

按参考文献所述（我还没太明白为什么作者说这里代进去的频率就是两个域的截止频率？），Buttworth滤波器的设计步骤如下：

1. 根据我们想要的数字滤波频率得到我们想要的模拟滤波器频率
2. 根据期望的模拟截止频率，将滤波器去归一化。
3. 进行双线性变换
4. 写成代码

Step1: 确定阶数和s域表达式

此处取一阶，N=1，所以s域表达式为：

$$H_a(s)=\frac{1}{1+s}\text{\hspace{1em}(3.10)}$$

Step2: 选取滤波器数字域的截止频率

假设为$f_d$

Step3: 计算对应s域的截止频率

根据式(3.8)可得：
$${\omega }_a=2\pi f_a=2f_s\tan (\frac{\pi f_d}{f_s})\text{\hspace{1em}(3.11)}$$

Step4: 去归一化

将式(3.10)中的s替代为：
$$s\to \frac{s}{{\omega }_a}$$
可得：
$$H_a(s)=\frac{1}{1+\frac{s}{{\omega }_a}}$$

Step5: 双线性变换

$$H(z)=\frac{Y(z)}{X(z)}=\frac{1}{1+\frac{2}{T{\omega }_a}\frac{1-z^{-1}}{1+z^{-1}}}=\frac{T{\omega }_a(1+z^{-1})}{T{\omega }_a(1+z^{-1})+2(1-z^{-1})}=\frac{T{\omega }_a(1+z^{-1})}{T{\omega }_a+2+(T{\omega }_a-2)z^{-1}}\text{\hspace{1em}(3.12)}$$
将式(3.11)代入上式，并设$\tan (\frac{\pi f_d}{f_s})=K$可得：
$$H(z)=\frac{Y(z)}{X(z)}=\frac{2K(1+z^{-1})}{2K+2+(2K-2)z^{-1}}=\frac{K(1+z^{-1})}{K+1+(K-1)z^{-1}}\text{\hspace{1em}(3.13)}$$
转换成差分方程：
$$\begin{gathered} \Rightarrow (K+1)y_k+(K-1)z^{-1}{y}_k=K(1+z^{-1})x_k \\ u\Rightarrow (K+1)y_k=K{x}_k+K{x}_{k-1}-(K-1)y_{k-1} \\ \Rightarrow y_k=\frac{K}{(K+1)}{x}_k+\frac{K}{(K+1)}{x}_{k-1}-\frac{(K-1)}{(K+1)}{y}_{k-1} \end{gathered}$$