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

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

数字系统的响应

数字滤波器是一种处理离散时间信号的系统。数字信号\(x[n]\)输入数字滤波器进行处理后生成输出信号\(y[n]\)。这个过程可以用卷积定理表示为：
\[y[n]=x[n]⊗h[n]=∑x[m]h[n-m] \tag{1}\] \(h[n]\)是数字系统在时域上的表示，称为冲激响应，也是系统输入信号为冲激序列时系统的输出结果。
对(1)做离散时间傅里叶变换，可以得到：
\[Y(e^{jω})=X(e^{jω})H(e^{jω}) \tag{2}\] \(H(e^{jω})\)是\(h[n]\)做离散时间傅里叶变换的结果，称为系统的频率响应。
> 可以发现(2)成立的条件中包含\(h[n]\)收敛的条件，如果\(h[n]\)不收敛，则不存在其离散时间傅里叶变换的结果\(H(e^{jω})\)。

数字系统的表示

除了上述两个表达式可以描述一个特定的数字系统外，一个特定的数字系统还可以使用如下两种方法进行表达：差分方程和传递函数。
系统可以表示为线性常系数差分方程：
\[∑_{k=0}^Na_ky[n-k]=∑_{k=0}^Mb_kx[n-k],a_0≠0,b_0≠0\] 对上式做Z变换，两边相比，得到系统的传递函数：
\[H(z)=\frac{Y(z)}{X(z)}=\frac{∑_{k=0}^Mb_kz^{-k}}{∑_{k=0}^Na_kz^{-k}}\] 但是由于z变换成立的条件是要给出收敛域RoC，因此在未给出传递函数的前提下，其对应的差分方程可能有两个： - 右边序列： \[∑_{k=0}^Na_ky[n-k]=∑_{k=0}^Mb_kx[n-k]\] 此时数字系统是因果系统。
- 左边序列: \[∑_{k=0}^Na_ky[n-k]=-∑_{k=-M}^{-1}b_kx[n+k]\] 此时数字系统是非因果系统。

传递函数的特征

与冲激响应

系统传递函数可以写作两个多项式的比：
\[H(z)=\frac{b_0}{a_0}\frac{Π_{k=1}^M(1-c_kz^{-1})}{Π_{k=1}^N(1-b_kz^{-1})}\] 其中\(M\)对应差分方程右侧除去\(x[n]\)外含有\(x\)项的个数，\(N\)对应差分方程左侧除去\(y[n]\)外含有\(y\)项的个数，\(c_k\)表示零点，\(d_k\)表示极点。
化简上式得到如下结构：
\[H(z)=∑_{r=0}^{M-N}B_rZ^{-r}+∑_{k=1}^N\frac{A_k}{1-d_kz^{-1}}\] 对其做Z反变换：
\[h[n]=∑_{r=0}^{M-N}B_rδ[n-r]+∑_{k=1}^NA_k(d_k)^nu[n]\] 对上述反变换结果进行讨论：
当\(N>0\)时，差分方程左边有除了\(y[n]\)的其他含\(y\)项\(y[n-k]\)，即此时系统方程当前输出由\(x\)项、系统以前的输出\(y[n-k]\)共同决定，系统含有反馈，此时系统方程为：
\[h[n]=∑_{r=0}^{M-N}B_rδ[n-r]+∑_{k=1}^NA_k(d_k)^nu[n]\] 由于\(u[n]\)是一个无限长度的序列，因此系统的冲激响应也是无限长度的，称这样的系统为无限冲激响应系统（IIR System）。
当\(N=0\)时，差分方程左边只含有\(y[n]\)，系统方程的当前输出只由输入\(x\)决定，系统不含有反馈，此时系统方程为：
\[h[n]=∑_{r=0}^{M}B_rδ[n-r]\] 由于\(δ[n-r]\)长度有限，整个系统的冲激响应的长度是有限的，称这样的系统为有限冲激响应系统（FIR System）。

	N	差分方程左侧	有无反馈	冲激响应长度
IIR	N>0	含有\(y[n-k]\)	有	无限
FIR	N=0	只有\(y[n]\)	无	有限

与频率响应

根据Z变换和离散时间傅里叶的定义式，可以推出：
\[H(e^{jω})=H(z)|_{z=exp(jω)}\] 因此系统的频率响应也可以写作多项式分式：
\[H(e^{jω})=\frac{b_0}{a_0}\frac{Π_{k=1}^M(1-c_ke^{-jω})}{Π_{k=1}^N(1-b_ke^{-jω})}\]

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

在\([0,2\pi]\)上生成若干个采样点，作为作图的\(x\)轴。
使用函数fft()可以生成指定序列的离散傅里叶变换结果。abs()函数能够给出指定序列在每个采样点上的绝对值，即幅度。angle()函数可以生成指定序列在每个采样点上的相位。
绘制其频率响应需要注意的是：
- 需要对\(h[n]\)进行周期延拓，方法是使用zeros(1,L)创建一个长度为L的全0向量用x=[x,zeros(1,L)]对其填充。
- 使用fft()函数求其FFT结果。abs()求结果的绝对值得到幅频响应，angle()求FFT结果的相频响应。
- 横轴应当为[0:(N+L)]*2*pi/(N+L+1)使得0-2π内每一格的长度为\(\frac{2π}{N+L+1}\) （因为0:(N+L)有N+L+1个数）
下面的示例程序展示了如何绘制\(h[n]=Sa(0.1(n-50)),0≤n≤100\)的幅频和相频响应曲线：

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{∑b_kz^{-k}}{∑a_kz^{-k}}\)时，可以使用[H,w]=freqz(b,a,n)返回其频率响应H和对应的π归一化的角频率横轴w，其中a是\(∑a_kz^{-k}\)由高次幂到低次幂排列时的系数向量，其中b是\(∑b_kz^{-k}\)由高次幂到低次幂排列时的系数向量，n为计算时所指定的N点FFT的点数，缺省值为512。
如下的例程中给出了如何绘制\(H(z)=\frac{1-0.5z^{-1}}{1-2z^{-1}+z^{-2}}\)的频率响应图：

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');