数字滤波器·数字系统的表示

数字系统的响应

数字滤波器是一种处理离散时间信号的系统。数字信号 $x[n]$ 输入数字滤波器进行处理后生成输出信号 $y[n]$ 。这个过程可以用卷积定理表示为：

$y[n] = x[n] \otimes h[n] = \sum x[m] h[n-m]$

$h[n]$ 是数字系统在时域上的表示，称为冲激响应，也是系统输入信号为冲激序列时系统的输出结果。对(1)做离散时间傅里叶变换，可以得到：

$Y(e^{j\omega}) = X(e^{j\omega}) H(e^{j\omega})$

$H(e^{j\omega})$ 是 $h[n]$ 做离散时间傅里叶变换的结果，称为系统的频率响应。

可以发现(2)成立的条件中包含 $h[n]$ 收敛的条件，如果 $h[n]$ 不收敛，则不存在其离散时间傅里叶变换的结果 $H(e^{j\omega})$ 。

数字系统的表示
除了上述两个表达式可以描述一个特定的数字系统外，一个特定的数字系统还可以使用如下两种方法进行表达：差分方程和传递函数。
系统可以表示为线性常系数差分方程：
$\sum_{k=0}^{N} a_{k} y[n-k] = \sum_{k=0}^{M} b_{k} x[n-k], a_{0} \neq 0, b_{0} \neq 0$
对上式做Z变换，两边相比，得到系统的传递函数：
$H(z) = \frac{Y(z)}{X(z)} = \frac{\sum_{k=0}^{M} b_{k} z^{-k}}{\sum_{k=0}^{N} a_{k} z^{-k}}$
但是由于z变换成立的条件是要给出收敛域RoC，因此在未给出传递函数的前提下，其对应的差分方程可能有两个：

右边序列： $\sum_{k=0}^{N} a_{k} y[n-k] = \sum_{k=0}^{M} b_{k} x[n-k]$ 此时数字系统是因果系统。

左边序列: $\sum_{k=0}^{N} a_{k} y[n-k] = - \sum_{k=-M}^{-1} b_{k} x[n+k]$ 此时数字系统是非因果系统。
传递函数的特征
与冲激响应
系统传递函数可以写作两个多项式的比： $H(z) = \frac{b_0}{a_0} \frac{\prod_{k=1}^{M}(1-c_kz^{-1})}{\prod_{k=1}^{N}(1-b_kz^{-1})}$ 其中 $M$ 对应差分方程右侧除去 $x[n]$ 外含有 $x$ 项的个数，$N$ 对应差分方程左侧除去 $y[n]$ 外含有 $y$ 项的个数，$c_k$ 表示零点，$d_k$ 表示极点。
化简上式得到如下结构： $H(z) = \sum_{r=0}^{M-N} B_r Z^{-r} + \sum_{k=1}^{N} \frac{A_k}{1-d_kz^{-1}}$ 对其做 Z 反变换： $h[n] = \sum_{r=0}^{M-N} B_r \delta[n-r] + \sum_{k=1}^{N} A_k (d_k)^n u[n]$ 对上述反变换结果进行讨论：
当 $N > 0$ 时，差分方程左边有除了 $y[n]$ 的其他含 $y$ 项 $y[n-k]$，即此时系统方程当前输出由 $x$ 项、系统以前的输出 $y[n-k]$ 共同决定，系统含有反馈，此时系统方程为： $h[n] = \sum_{r=0}^{M-N} B_r \delta[n-r] + \sum_{k=1}^{N} A_k (d_k)^n u[n]$ 由于 $u[n]$ 是一个无限长度的序列，因此系统的冲激响应也是无限长度的，称这样的系统为无限冲激响应系统（IIR System）。
当 $N = 0$ 时，差分方程左边只含有 $y[n]$，系统方程的当前输出只由输入 $x$ 决定，系统不含有反馈，此时系统方程为： $h[n] = \sum_{r=0}^{M} B_r \delta[n-r]$ 由于 $\delta[n-r]$ 长度有限，整个系统的冲激响应的长度是有限的，称这样的系统为有限冲激响应系统（FIR System）。
| | N | 差分方程左侧 | 有无反馈 | 冲激响应长度 |
| —- | —- | —- | —- | —- |
| IIR | N>0 | 含有 $y[n-k]$ | 有 | 无限 |
| FIR | N=0 | 只有 $y[n]$ | 无 | 有限 |
与频率响应
根据Z变换和离散时间傅里叶的定义式，可以推出： $H(e^{j\omega}) = H(z)|_{z=exp(j\omega)}$ 因此系统的频率响应也可以写作多项式分式： $H(e^{j\omega}) = \frac{b_0}{a_0} \frac{\Pi_{k=1}^{M}(1-c_ke^{-j\omega})}{\Pi_{k=1}^{N}(1-b_ke^{-j\omega})}$

使用MATLAB® 绘制系统的频率响应

在 [0, 2π] 上生成若干个采样点，作为作图的 x 轴。

使用函数 fft() 可以生成指定序列的离散傅里叶变换结果。 abs() 函数能够给出指定序列在每个采样点上的绝对值，即幅度。 angle() 函数可以生成指定序列在每个采样点上的相位。
绘制其频率响应需要注总的是：

·需要对$h[n]$进行周期延拓，方法是使用 zeros(1,L)创建一个长度为L的全0向量用 x=[x,zeros(1,L)] 对

其填充。

·使用 fft()函数求其FFT结果。 abs()求结果的绝对值得到幅频响应，angle()求FFT结果的相频响

应。

· 横轴应当为 [0:(N+L)]2pi/(N+L+1) 使得O-2π内每一格的长度为$\frac{2\pi}{N+L+1}$ (因为 $\theta:(\mathbb{N}+L)$ 有N+L+1个

数)

下面的示例程序展示了如何绘制$h[n]=Sa(0.1(n-50)),0\leq n\leq100$的幅频和相频响应曲线：
```shell
N=100; %100点FFT
n=0:N;
x=sinc(0.1(n-50));
x=[x,zeros(1,1000)]; %对x周期延拓到长度为1000
X=fft(x); %求其FFT/DFS结果（Amplitude Response）
subplot(2,1,1);
w=[0:(N+1000)]2*pi/(N+1001);%横轴，每一格表示2π/(N+1001)
plot(w./pi,abs(X)); %取绝对值求幅频响应（Magnitude Response）并作π归一化
title(‘Magnitude Response’);
axis ([0 2 0 12]);
xlabel(‘\omega/\pi’);
subplot(2,1,2);
plot(w./pi,angle(X)); % angle(X)表示求其相频响应
title(‘Phase Response’);
axis ([0 2 -4 4]);
xlabel(‘\omega/\pi’);

```

给定系统Z域下的传递函数$H(z)=\frac{\sum b_kz^{-k}}{\sum a_kz^{-k}}$时，可以使用[H,w]=freqz(b,a,n)返回其频率响应 н和对应的π归一化的角频率横轴 w,其中 a 是$\sum a_kz^{-k}$由高次幂到低次幂排列时的系数向量，其中 b 是$\sum b_kz^{-k}$由高次幂到低次幂排列时的系数向量，n 为计算时所指定的N点FT的点数，缺省值为512。
如下的例程中给出了如何绘制$H(z)=\frac{1-0.5z^{-1}}{1-2z^{-1}+z^{-2}}$的频率响应图：
```shell
b=[1,-0.5];
a=[1,-2,1];
[H,w]=freqz(b,a);
plot(w,abs(X)); %取绝对值求幅频响应（Magnitude Response）
title(‘Magnitude Response’);
plot(w,angle(X)); % angle(X)表示求其相频响应
title(‘Phase Response’);
```