滤波器设计实验

姓名：周镇峰
学号：201830260422
班级：2018电子科学与技术卓越班

数字信号处理实验三

实验目的

• 利用窗函数法和等波纹滤波器法设计 FIR 滤波器
• 利用脉冲响应不变法和双线性变换法设计 IIR 滤波器
• 利用 STFT 进行参数分析

实验内容

1. 利用 Matlab 编程，分别用窗函数法和等波纹滤波器法设计两种 FIR 数字带通滤波器，指标要求如下：
   • 通带边缘频率 $\Omega_{P1}=0.45\pi$，$\Omega_{P2}=0.65\pi$，通带峰值起伏：$\alpha_P \leq 1$ dB
   • 阻带边缘频率 $\Omega_{S1}=0.3\pi$，$\Omega_{S1}=0.75\pi$，最小阻带衰减 $\alpha_S \geq 40$ dB
2. 利用 Matlab 编程，用脉冲响应不变法和双线性变换法设计一个数字带通滤波器，指标要求如下：
   • 通带边缘频率 $\Omega_{P1}=0.45\pi$，$\Omega_{P2}=0.65\pi$，通带峰值起伏：$\alpha_P \leq 1$ dB
   • 阻带边缘频率 $\Omega_{S1}=0.3\pi$，$\Omega_{S1}=0.8\pi$，最小阻带衰减 $\alpha_S \geq 40$ dB
3. 线性调频信号LFM $x[nT]=\exp (j2\pi f_c nT + j \pi \mu (nT)^2)$，要求
   1. 通过短时傅里叶变换，求给定起始频率 $f_c$ 和调频率 $\mu$ 的 LFM 信号的时频图
   2. 查找其他参数估计方法，通过时频图重新估计 LFM 信号的起始频率和调频率

实验过程

设计 FIR 滤波器

窗函数法

  窗函数法就是利用一个窗函数+时移，将无穷、非因果的理想滤波器转换为物理可实现的滤波器。加窗的过程如下图所示：

  常见的窗函数的特性见下表：

  根据题目要求，最小阻带衰减 $\alpha_S \geq 40$ dB，因此可以选择汉宁窗。

  过渡带 $\Delta \Omega_1 = \Omega_{P1}-\Omega_{S1}=0.15\pi$，$\Delta \Omega_2 = \Omega_{S2}-\Omega_{P2}=0.1\pi$，取两者中较小的，代入上表中 $\Delta \Omega_2=3.11\pi/M$，可得 $M\approx 32$

  利用 hamming(M) 函数生成 $M$ 长的汉明窗，并利用 fir1() 完成基于窗函数的 FIR 滤波器设计

wp=[0.45 0.65]; %通带
ws=[0.3 0.75]; %阻带
wc=(wp+ws)/2; %截止频率
as=40;
w=min(abs(wp-ws)); %过渡带
M=ceil(3.11/w); %长度
win=hamming(M+1);
boxb=fir1(M,wc,win);
freqz(boxb,1,512);

等波纹滤波器法

  等波纹最佳逼近法是一种优化设计法，即最大误差最小化准则，它克服了窗函数设计法和频率采样法的缺点，使最大误差（即波纹的峰值）最小化，并在整个逼近频段上均匀分布。用等波纹最佳逼近法设计的FIR数字滤波器的幅频响应在通带和阻带都是等波纹的，而且可以分别控制通带和阻带波纹幅度，这就是等波纹的含义。最佳逼近是指在滤波器长度给定的条件下，使加权误差波纹幅度最小化。与窗函数设计法和频率采样法比较，由于这种设计法使滤波器的最大逼近误差均匀分布，所以设计的滤波器性能价格比最高。阶数相同时，这种设计法使滤波器的最大逼近误差最小，即通带最大衰减最小，阻带最小衰减最大；指标相同时，这种设计法使滤波器阶数最低。

  等波纹最佳逼近法的设计思想:用 $H_d(\omega)$ 表示希望逼近的幅度特性函数，要求设计线性相位FIR数字滤波器时,$H_d(\omega)$ 必须满足线性相位约束条件。用 $H(\omega)$ 表示实际设计的滤波器的幅度特性函数。定义加权误差函数$\varepsilon(\omega)$为：

  式中，$W(\omega)$ 为幅度误差加权函数，用来控制不同频带（一般指通带和阻带）的幅度逼近精度。等波纹最佳逼近法的设计在于找到滤波器的系数向量 $h[n]$，使得在通带和阻带内的最大绝对值幅度误差 $|\varepsilon(\omega)|$为最小，这也就是最大误差最小化问题。

  在 MATLAB 中，调用MATLAB信号处理工具箱函数remezord来计算等波纹滤波器阶数N和加权函数W(ω)，调用函数remez可进行等波纹滤波器的设计，直接求出滤波器系数。函数remezord中的数组fedge为通带和阻带边界频率，数组mval是两个边界处的幅值，而数组dev是通带和阻带的波动，fs是采样频率单位为Hz。

[n,fpts,mag,wt]=remezord([0.3 0.45 0.65 0.75],[0 1 0],[0.01 0.1087 0.01]);
%用remezord函数估算出remez函数要用到的阶n、归一化频带边缘矢量fpts、频带内幅值响应矢量mag及加权矢量w，
%使remez函数设计出的滤波器满足f、a及dev指定的性能要求。
h2=remez(n,fpts,mag,wt);%设计出等波纹滤波器
[hh2,w2]=freqz(h2,1,256);
figure(2)
subplot(2,1,1)
plot(w2/pi,20*log10(abs(hh2)))
grid
xlabel('归一化频率w');ylabel('幅度/db');
subplot(2,1,2)
plot(w2/pi,angle(hh2))
grid
xlabel('归一化频率w');ylabel('相位/rad');

设计 IIR 滤波器

  总的来说，都是利用模拟滤波器来设计数字滤波器，根据转换的方法不同，分为：

• 冲激响应不变法
• 双线性变换法

冲激响应不变法

  变换原理：数字滤波器的脉冲响应为 $h[n]$，模拟滤波器的冲激响应为 $h_a(t)$，冲激响应不变就是指这俩在 $t=nT$ 时相同，即：

  映射关系如下图（多对一）：

  由于数字滤波器相当于对模拟滤波器在时域上采样，频域上进行周期延拓，所以只适合设计低通和带通滤波器（避免混叠）。这种方法的优点是时域逼近良好，且相位保持线性关系。

  模拟到数字的转换过程如下：

• $H_a(s)=\sum_{k=1}^N \dfrac{A_k}{s-p_k}$
• $h_a(t)=\sum_{k=1}^N A_k e^{p_k t}u(t)$
• $h(n)=h_a(nT)=\sum_{k=1}^N A_k (e^{p_kT})^n u(n)$
• $H(z)=\sum_{n=-\infty}^\infty h(n)z^{-n}$ $=\sum_{n=0}^{\infty}\sum_{k=1}^N A_k \left( e^{p_kT}\right)^n z^{-n}$ $=\sum_{k=1}^N \frac{A_k}{1-e^{p_k T}z^{-1}}$

  最终我们得到：

  综上，设计过程是这样的：

1. 选择 $T$，然后求出对应的模拟频率： $$ \Omega = \frac{\omega}{T} $$
2. 根据 $\Omega_p$，$\Omega_s$，$R_p$，$A_s$ 设计模拟滤波器
3. 将 $H_a(s)$ 展开为部分分式的形式： $$ H_a(s)=\sum_{k=1}^{N} \frac{R_k}{s-p_k} $$
4. 将模拟转化为数字： $$ H(z)=\sum_{k=1}^N \frac{R_k}{1-e^{p_k T}z^{-1}} $$

  我们可以用 MATLAB 来辅助我们设计模拟滤波器。下面代码中的 buttord 和 butter 就是设计巴特沃斯滤波器。

T=1; %采样周期
fs=1/T; %采样频率
wp=[0.45*pi 0.65*pi]*fs; %通带
Ap=1; %通带峰值起伏dB
ws=[0.3*pi 0.8*pi]*fs; %阻带
As=40; %阻带峰值起伏dB
[N, wc]=buttord(wp,ws,Ap,As,'s'); %得到 buttonworth 滤波器的阶数N和截止频率wc
[B, A]=butter(N,wc,'bandpass','s'); %设计 buttonworth 滤波器，得到H(s)=B/A
w=linspace(0,pi,400*pi);
[D,C]=impinvar(B,A,fs); %模拟转数字
Hz=freqz(D,C,w);
plot(w/pi,abs(Hz));
grid on;
title('巴特沃斯带通滤波器');
xlabel('Frequency/Hz');
ylabel('Magnitude');

  我们也可以使用其他滤波器，比如切比雪夫I型滤波器：

T=1; %采样周期
fs=1/T; %采样频率
wp=[0.45*pi 0.65*pi]*fs; %通带
Ap=1; %通带峰值起伏dB
ws=[0.3*pi 0.8*pi]*fs; %阻带
As=40; %阻带峰值起伏dB
[N, wc]=cheb1ord(wp,ws,Ap,As,'s'); %得到 buttonworth 滤波器的阶数N和截止频率wc
[B, A]=cheby1(N,Ap, wc,'bandpass','s'); %设计 buttonworth 滤波器，得到H(s)=B/A
w=linspace(0,pi,400*pi);
[D,C]=impinvar(B,A,fs); %模拟转数字%冲激响应不变
Hz=freqz(D,C,w);
plot(w/pi,abs(Hz));
grid on;
title('切比雪夫Ⅰ型带通滤波器');
xlabel('Frequency/Hz');
ylabel('Magnitude');

双线性变换法

  双线性变换就是用了不同的 $s\rightarrow z$ 转换方法：

  它的映射关系如下（一对一）：

  这种变换的优点就是：s平面与z平面是单值变换，所以避免了混叠。缺点就是：除了零频率附近，模拟频率 $\Omega$ 与数字频率 $\omega$ 之间非线性，因此要求幅度响应应该为分段常数型，否则会产生畸变。

  公式的推导过程如下：

  我们先将模拟频率 $\Omega=[-\infty,+\infty]$ 压缩到 $\Omega_1=[-\pi/T，\pi/T]$

  然后再将 $\Omega_1$ 映射到 $z$（就是冲激响应不变法中的过程）：

  将两个结合起来，有：

  推导完毕。从上面推导中，可以注意到模拟频率 $\Omega$ 与数字频率 $\omega=\Omega_1 T$ 的对应关系不是线性的，而我们希望它越接近线性越好，于是我们进行一个拉伸操作：

  若要求在低频处有比较明确的（线性）对应关系，那么：

  从而双线性变换可以进一步写为：

  这正是课本中介绍的双线性变换，也是下面程序中用到的双线性变换。

  若要求在某个频率 $\Omega_c\leftrightarrow \omega_c$ 处有比较明确的（线性）对应关系，那么：

  式中的 $\Omega_c, \omega_c$ 是自己设定的。

  MATLAB 实现的代码如下：

T=1; %采样周期
fs=1/T; %采样频率
wp=[0.45*pi 0.65*pi]; %通带
Ap=1; %通带峰值起伏dB
ws=[0.3*pi 0.8*pi]; %阻带
As=40; %阻带峰值起伏dB

wp=(2/T)*tan(wp/2); %预扭曲%数字频率转模拟频率
ws=(2/T)*tan(ws/2);

[N, wc]=buttord(wp,ws,Ap,As,'s'); %得到 buttonworth 滤波器的阶数N和截止频率wc
[B, A]=butter(N,wc,'bandpass','s'); %设计 buttonworth 滤波器，得到H(s)=B/A
w=linspace(0,pi,400*pi);
[D,C]=bilinear(B,A,fs); %模拟转数字%双线性
Hz=freqz(D,C,w);
plot(w/pi,abs(Hz));
grid on;
title('巴特沃斯带通滤波器');
xlabel('Frequency/Hz');
ylabel('Magnitude');

利用 STFT 对 LFM 信号进行参数估计

可以利用 chirp() 函数生成 LFM 信号：

fs=1e3; %采样率
T=2; %脉冲宽度
N=fs*T; %采样点数
t=linspace(0, T, N);
f0=100; %起始频率
f1=400; %终止频率
y=chirp(t, f0, T, f1);
plot(t,y);

利用 spectrogram() 计算短时傅里叶变换：

spectrogram(y,128,100,128,1e3,'yaxis');

取图像中每一列的最大值，作出图像：

[~,f,t,p] = spectrogram(y,128,100,128,1e3,'yaxis');
[fridge,~,lr] = tfridge(p,f);
plot3(t,fridge,abs(p(lr)),'LineWidth',4)
t(1),fridge(1) %第一个点的坐标
t(end),fridge(end) %最后一个点的坐标

图中，首尾的坐标为：$(0.064, 109.375)$，$(1.912,390.625)$，可以估算出斜率为：$\frac{390.625-109.375}{1.912-0.064}=152.1916$，实际的斜率为 $\frac{400-100}{2-0}=150$，误差为 $1.46\%$

通过上面的实验，STFT 法简单易用，但存在两个困难：1. 窗函数的选择；2. 窗长度的选择。根据 Heisenberg 测不准原则，信号的时/频分辨力不能同时提高。

实验心得

通过 matlab 实验，加深了对几种滤波器设计方法的了解。同时，粗略地了解了利用 STFT 进行参数估计的方法。

参考

• 课本
• 用MATLAB结合窗函数法设计数字带通FIR滤波器
• 冲激响应不变法与双线性变换
• 用冲激响应不变法设计IIR数字滤波器