滤波器设计实验

\(\begin{align*} \newcommand{\dif}{\mathop{}\!\mathrm{d}} \newcommand{\belowarrow}[1]{\mathop{#1}\limits_{\uparrow}} \newcommand{\bd}{\boldsymbol} \newcommand{\L}{\mathscr{L}} \newcommand{\xleftrightarrow}[1]{\stackrel{#1}{\longleftrightarrow}} \end{align*}\)

姓名:周镇峰
学号: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 滤波器

窗函数法

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

FIR 窗函数法

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

窗类型 主瓣宽度 过渡带宽度 相对旁瓣水平 最小阻带衰减
矩形窗 4π/(2M+1) 0.92π/M -13.3dB -20.9dB
汉宁窗 8π/(2M+1) 3.11π/M -31.5dB -43.9dB
汉明窗 8π/(2M+1) 3.32π/M -42.7dB -54.5dB
布莱克曼窗 12π/(2M+1) 5.56π/M -58.1dB -75.3dB

  根据题目要求,最小阻带衰减 $\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

等波纹滤波器法

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

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

\[\varepsilon(\omega)=W(\omega)[H_d(\omega)-H(\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');

等波纹法设计FIR

设计 IIR 滤波器

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

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

冲激响应不变法

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

\[h[n]=h_a(t) \Big\vert_{t=nT}\\ H(z)\Big\vert_{z=e^{sT}} = \hat{H}_a (s)=\frac{1}{T}\sum_{k=-\infty}^{+\infty} H_a \left( s-j\frac{2\pi}{T}k \right)\]

  映射关系如下图(多对一):

冲激响应不变法的映射关系

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

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

\[H_a(s) \rightarrow h_a(t) \rightarrow h_a(nT) \rightarrow h(n) \rightarrow H(z)\]
  • $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}}$

  最终我们得到:

\[H_a(s)=\sum_{k=1}^N \frac{A_k}{s-p_k} \rightarrow H(z)=\sum_{k=1}^N \frac{A_k}{1-e^{s_kT}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 来辅助我们设计模拟滤波器。下面代码中的 buttordbutter 就是设计巴特沃斯滤波器。

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

dsp巴特沃斯

  我们也可以使用其他滤波器,比如切比雪夫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');

dsp切比雪夫1型

双线性变换法

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

\[s=\frac{1-z^{-1}}{1+z^{-1}}\\ z=\frac{1+s}{1-s}\]

  它的映射关系如下(一对一):

双线性变换法的映射关系

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

  公式的推导过程如下:

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

\[\Omega=\tan \frac{\Omega_1 T}{2}\]

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

\[z=e^{s_1 T}\]

  将两个结合起来,有:

\[\begin{aligned} \Omega&=\tan \frac{\Omega_1 T}{2}\\ &=\frac{\sin\frac{\Omega_1 T}{2}}{\cos\frac{\Omega_1 T}{2}}\\ &=\frac{e^{j\frac{\Omega_1 T}{2}}-e^{-j\frac{\Omega_1 T}{2}}}{2j}\bigg/ \frac{e^{j\frac{\Omega_1 T}{2}}+e^{-j\frac{\Omega_1 T}{2}}}{2} \end{aligned}\] \[\begin{aligned} s&=j\Omega\\ &=\frac{e^{j\frac{\Omega_1 T}{2}}-e^{-j\frac{\Omega_1 T}{2}}}{e^{j\frac{\Omega_1 T}{2}}+e^{-j\frac{\Omega_1 T}{2}}}\\ &=\frac{1-e^{-s_1 T}}{1+e^{-s_1T}}\\ &=\frac{1-z^{-1}}{1+z^{-1}} \end{aligned}\]

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

\[\Omega = c \cdot \tan \frac{\omega}{2}\\ s = c \frac{1-z^{-1}}{1+z^{-1}}\\ z=\frac{c+s}{c-s}\]

  若要求在低频处有比较明确的(线性)对应关系,那么:

\[\Omega_1 \approx \Omega = c \cdot \tan \frac{\Omega_1 T}{2} \approx c \cdot \frac{\Omega_1T}{2}\\ \Rightarrow c = \frac{2}{T}\]

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

\[s = \frac{2}{T} \left(\frac{1-z^{-1}}{1+z^{-1}}\right)\]

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

  若要求在某个频率 $\Omega_c\leftrightarrow \omega_c$ 处有比较明确的(线性)对应关系,那么:

\[\Omega_c = c \cdot \tan \frac{\omega_c}{2}\\ \Rightarrow c = \Omega_c \cot \frac{\omega_c}{2}\]

  式中的 $\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');

dsp巴特沃斯_双线性.jpg

利用 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);

chirp信号

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

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

stft

取图像中每一列的最大值,作出图像:

[~,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) %最后一个点的坐标

stft最大值

图中,首尾的坐标为:$(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 进行参数估计的方法。

参考