DSP第三次上机作业

初夏十九

于 2022-12-02 11:44:34 发布

阅读量278

点赞数 2

分类专栏： DSP上机作业文章标签： matlab

本文链接：https://blog.csdn.net/weixin_54586334/article/details/128145476

版权

DSP上机作业专栏收录该内容

3 篇文章 1 订阅

订阅专栏

DSP第三次上机作业

一、实验目的

学习使用MATLAB设计数字滤波器，掌握IIR和FIR数字滤波器的设计方法。

二、Task 1

2.1 题目要求

（1）用双线性变换法设计一个数字低通Butterworth滤波器，以满足：
通带截止频率： $\omega_p=0.2\pi$ ,通带最大衰减： $R_p=1dB$
阻带截止频率： $\omega_s=0.3\pi$ ,阻带最小衰减： $R_s=15dB$
（2）绘制其频率响应的幅频特性与相频特性图
（3）通过一个含有低频与高频成份的混合信号进行滤波验证：
$x(t)=\sin{(2\pi f_1t)}+\sin{(2\pi f_2t)}, f_1=50Hz,f_2=250Hz$

2.2 设计思路

双极性变换法设计数字滤波器步骤为：
（1）确定参数 $T$
（2）将数字滤波器的边界频率 $\omega_p$ 、 $\omega_s$ 转换为模拟滤波器的边界频率 $\Omega_p$ 、 $\Omega_s$ ，即预畸变
（3）按照模拟滤波器技术指标设计模拟滤波器 $H_a(s)$
（4）用双极性变换法将模拟滤波器变换为数字滤波器
由于 $T$ 的选择与最终设计出的数字滤波器无关，为方便计算，我们令 $T = 2$ 。程序采用的方法是先设计巴特沃斯模拟低通滤波器，然后利用双线性变换函数，将模拟滤波器变换为数字滤波器。在第三问中，采样频率设置为1000Hz，采样点设置100个，由此绘制时域波形和频域波形。

2.3 程序源码

%% 双线性变换法设计数字滤波器
clear;clc;
Wp=0.2*pi;      %通带截止频率
Ws=0.3*pi;      %阻带截止频率
Rp=1;           %通带最大衰减
Rs=15;          %阻带最小衰减
T=2;            %令T=2，T的选择与最终设计出的数字滤波器无关

%预畸变
Omega_p=2/T*tan(Wp/2);       
Omega_s=2/T*tan(Ws/2);

%计算巴特沃斯模拟低通滤波器的阶数和3dB截止频率
[N,Omega_c]=buttord(Omega_p,Omega_s,Rp,Rs,'s'); 
[bs,as]=butter(N,Omega_c,'s');  %计算巴特沃斯模拟滤波器系统函数中分子和分母多项式系数向量B和A
[bz,az]=bilinear(bs,as,1/T);    %用双线性变换法将模拟滤波器变换为数字滤波器

%% 绘制幅频特性与相频特性曲线
W=linspace(0,pi,2000);          %在[0,pi]区间内选取2000个采样点
H=freqz(bz,az,W);               %Z域频率响应，为复数形式

subplot(2,3,1);
plot(W/pi,20*log10(abs(H)));    %绘制系统幅频特性曲线（单位：dB）
title('幅频特性');xlabel('\omega/\pi');ylabel('|H(z)|/dB')
subplot(2,3,4);
plot(W/pi,angle(H)*180/pi);     %绘制系统相频特性曲线（单位：Deg）
title('相频特性');xlabel('\omega/\pi');ylabel('\phi[H(z)]/Deg')

%% 通过混合信号验证滤波功能
fs=1000;  %采样频率
N=100;    %采样点100个
n=0:N-1;
t=n/fs;   %t=nT
f=fs/N*n;

x=sin(2*pi*50*t)+sin(2*pi*250*t);  
subplot(2,3,2);         %绘制输入信号时域曲线
plot(t,x);title('输入信号时域波形');xlabel('t/ms');ylabel('x')

X=fft(x,N);Xm=abs(X);
subplot(2,3,5);         %绘制输入信号频域曲线
plot(f(1:N/2),Xm(1:N/2)); title('输入信号频域波形');xlabel('f/Hz');ylabel('X(k)')

y=filter(bz,az,x);
subplot(2,3,3);         %绘制输出信号时域曲线
plot(t,y);title('输出信号时域波形');xlabel('t/ms');ylabel('y')

Y=fft(y,N);Ym=abs(Y);
subplot(2,3,6);         %绘制输出信号频域曲线
plot(f(1:N/2),Ym(1:N/2)); title('输出信号频域波形');xlabel('f/Hz');ylabel('Y(k)')

2.4 程序结果

数字滤波器的幅频特性和相频特性：
Task1结果图1
输入、输出信号的时域、频域波形：
Task1结果图2

2.5 结果分析

从幅频特性图中可以看出，所设计的数字滤波器满足题目中给出的各项指标要求。其中通带最大衰减约为-0.62dB，小于要求的1dB，阻带最小衰减约为-15.13dB，大于要求的15dB。从相频特性图中可以看出，滤波器的相频响应为非线性曲线.由于angle函数的取值是 $-\pi$ 到 $\pi$ ，故当相移下降到 $- 180 °$ 时，会跳变至 $+ 180 °$ 后再开始下降，实际上相移会随着 $\omega$ 的不断增大，从 $0 °$ 一直下降到 $- 360 °$ 。
在利用混合信号验证滤波器功能时，输入混合信号为两个频率分别为50kHz和250kHz的正弦信号之和，从频谱图中可以清楚看到50kHz和250kHz两个幅度尖峰。因为程序设置的采样频率为1000Hz，所以根据公式 $\omega=\Omega T$ ，经过离散化采样后信号的数字频率分别为 $0.1\pi$ 和 $0.5\pi$ ，一个位于滤波器的通带内，一个位于滤波器的阻带内。理论上，信号通过数字滤波器后，位于阻带内的250kHz频率分量将被滤除掉，输出信号只保留50kHz的频率分量。从时域波形图可以看出，输出信号为50kHz的单音正弦信号，频谱图中只有50kHz一个幅度尖峰，与理论相符合。程序结果说明作者设计的滤波器实现了低通滤波的功能，效果良好。

三、Task 2

3.1 题目要求

（课本习题7.13）希望对输入低频模拟信号做数字低通滤波处理，以提取所需要的信号。设系统的采样频率为 $50 k H z$ ，要求通带截止频率为 $10 k H z$ ，阻带截止频率为 $25 k H z$ ，阻带最小衰减为 $60 d B$ 。
（1）用窗函数设计法设计FIR数字低通滤波器，选择合适的窗函数及滤波器长度N，求出单位脉冲响应 $h (n)$ 。
（2）画出该FIR数字低通滤波器的幅频响应特性曲线和相频响应特性曲线。

3.2 设计思路

窗函数设计法设计FIR数字滤波器的步骤为：
（1）根据所需设计的数字滤波器类型，确定线性相位数字滤波器的类型（Ⅰ、Ⅱ、Ⅲ、Ⅳ型）
（2）根据滤波器阻带衰减 $\alpha_s$ ，选择合适的窗函数 $\omega(n)$
（3）确定理想数字滤波器的频率响应函数 $H_d(e^{j\omega})=H_d(\omega)\ e^{j\theta_d(\omega)}$ ，其中 $H_d(\omega)$ 为幅度特性函数， $\theta_d(\omega)$ 为相位特性函数
（4）计算理想数字滤波器的单位脉冲响应 $h_d(n)$
（5）加窗得到设计结果 $h (n)$ ,即 $h(n)=h_d(n)\omega(n)$ 题目中要求阻带衰减 $\alpha_s$ =60dB，故应选择布莱克曼窗。程序采用I型滤波器设计低通滤波器，由布莱克曼窗过渡带宽度 $∆B=\frac{11\pi}{N}\le\frac{(25-10)\times2\pi\times{10}^3}{50\times{10}^3}$ 解得 $N\geq18.33$ Ⅰ型滤波器长度为奇数，故取N=19。

3.3 程序源码

clc,clear;
fs=50000;               %采样频率50kHz
Omega_p=2*pi*10000;     %通带截止频率10kHz
Omega_s=2*pi*25000;     %阻带截止频率25kHz
Wp=Omega_p/fs;
Ws=Omega_s/fs;
Wc=(Wp+Ws)/2/pi;        %计算对pi归一化的6dB数字截止频率

N=19;
n=0:N-1;
h=fir1(N-1,Wc,blackman(N));     %布莱克曼窗
subplot(3,1,1);                 %绘制数字滤波器的单位脉冲响应
stem(n,h,'.');xlabel('n');ylabel('h(n)')

W=linspace(0,pi,512);           %在[0,pi]区间内选取512个采样点
H=freqz(h,1,W);                 %Z域频率响应
subplot(3,1,2);
plot(W/pi,20*log10(abs(H)));    %绘制系统幅频特性曲线（单位：dB）
title('幅频特性');xlabel('\omega/\pi');ylabel('|H(z)|/dB')
subplot(3,1,3);
plot(W/pi,angle(H)*180/pi);     %绘制系统相频特性曲线（单位：Deg）
title('相频特性');xlabel('\omega/\pi');ylabel('\phi[H(z)]/Deg')

3.4 程序结果

数字滤波器的单位脉冲响应：
Task2结果图1
数字滤波器的幅频特性和相频特性：
Task2结果图2

3.5 结果分析

题目中给定采样频率为50kHz，因此根据 $\omega =\Omega T$ ，阻带截止频率对应的数字频率为 $\omega=1\pi$ ，从幅频特性图可以看出，数字滤波器阻带截止频率对应的衰减值为-65.06dB，小于要求的-60dB，满足设计要求。从相频特性图中可以看出，所设计的数字滤波器为线性相位滤波器，相移从 $0 °$ 一直下降到 $- 1620 °$ 。

四、Task 3

4.1 题目要求

（课本习题7.30）利用MATLAB工具箱函数并编写程序，采用频率采样法设计一个严格线性相位FIR数字低通滤波器，要求通带截止频率 $\omega_p=\pi/3rad$ ，过渡带宽度 $∆_B\le\pi/16rad$ ，阻带最小衰减 $\alpha_s=40dB$ ，允许在过渡带上设置采样点。画出所设计滤波器的单位脉冲响应 $h (n)$ 的波形及幅频响应特性曲线。

4.2 设计思路

频率采样法设计FIR数字滤波器的步骤为：
（1）根据阻带最小衰减 $\alpha_s$ ，确定过渡带采样点的个数m
（2）根据过渡带宽度 $\Delta B$ 的要求，估算滤波器的长度N
（3）构造希望逼近滤波器的频率响应函数
（4）对 $H_d(e^{j\omega})$ 进行频域等间隔N点采样，得 $H(k)={H_d(e^{j\omega})|}_{\omega=\frac{2\pi}{N}k}=H_ke^{j\theta_k},\ \ 0 \le k \le1$ 其中 $H_k=H_d(\omega)|_{\omega=\frac{2\pi}{N}k}$ ， $\theta_k=\theta(\omega)|_{\omega=\frac{2\pi}{N}k}$
（5）对H（k）进行N点IDFT，得所设计FIR数字滤波器的单位脉冲响应 $h (n)$
题目中要求 $\alpha_s=40dB$ ，故设置1个过渡带采样点，由 $N\geq\frac{(m+1)2π}{\Delta B}$ 计算得到滤波器的阶数 $N\geq64$ ，程序采用Ⅰ型滤波器设计低通滤波器，故取N=65。希望逼近滤波器选择理想低通滤波器，即 $H_d(e^{j\omega})=H_d(\omega)\ e^{-j\omega(N-1)/2}$ 。对 $H_d(e^{j\omega})$ 的频域等间隔N点采样是在 $[0,2\pi]$ 区间上，对应归一化的 $\omega$ 取值范围为[0,2],由于N为奇数，故 $H_k$ 偶对称于 $\omega=\pi$ 。通带截止频率 $\omega_p=\pi/3rad$ ，满足 $\frac{21}{65}rad\le\omega_p\le\frac{22}{65}rad$ ，故得 $\theta_k=-\omega\frac{N-1}{2}=-k\frac{2\pi}{N}\frac{N-1}{2} =-\frac{64}{65}\pi k$ $H_k=\left\{\begin{aligned} 1,\ \ 0\le k\le11;56\le k\le65 \\ H_b,\ \ k=12;k=55 \\ 0,\ \ 13\le k\le54 \end{aligned} \right.$ 其中 $H_b$ 为过渡点采样值，在 $k = 12$ ， $k = 55$ 处安排过渡点采样。

4.3 程序源码

clc;clear;
N=65;       %滤波器长度65
n=0:64;
N1=11;      
N2=N-N1-(N1-1)-2; 
H_b=0.38;   %过渡带采样值
Theta_k=-pi*64/65*n;
H_k=[ones(1,N1),H_b,zeros(1,N2),H_b,ones(1,N1-1)];
subplot(3,1,1);                 %绘制频率采样后的幅度特性函数
stem(n,H_k,'.');xlabel('\omega/\pi');ylabel('Hk')

H=H_k.*exp(1i*Theta_k);
h=real(ifft(H));
subplot(3,1,2);                 %绘制数字滤波器的单位脉冲响应
stem(n,h,'.');title('单位脉冲响应');xlabel('n');ylabel('h(n)')

W=linspace(0,pi,512);           %在[0,pi]区间内选取512个采样点
H=freqz(h,1,W);                 
subplot(3,1,3);
plot(W/pi,20*log10(abs(H)));    %绘制系统幅频特性曲线（单位：dB）
title('幅频特性');xlabel('\omega/\pi');ylabel('|H(z)|/dB')