用matlab设计高通滤波器,雪比切夫、fir两种方法 课程设计HPF

课 程 设 计
20011 年 7月 1日 设计题目 学
号
专业班级 指导教师 学生姓名 张腾达 吴晔 陈丽娟 杨蕾
通信电子电路课程设计 ——数字滤波器的设计 张静 光信息08-3班
实验组员 张静 胡磊 艾永春 赵亚龙
王宏道 胡进娟 马丽婷
设计要求：
某系统接收端接收到的信号为
y=cos(2π*60t)+1.2cos(2π*140t)+2sin(2π*220t)
+1.5sin(2π*300t)
(A) 发现此信号夹杂了一个正弦噪声noise=1.5sin(2π*300t)，请设计一个低通滤波器将此噪声滤除，从而恢复原信号。

(B) 发现此信号夹杂了一个正弦噪声noise= cos(2π*60t)
+1.5sin(2π*300t) ，请设计一个带通滤波器将此噪声滤除，从而恢复原信号。

(C) 发现此信号夹杂了一个正弦噪声noise=
1.2cos(2π*140t)+2sin(2π*220t)，请设计一个带阻滤波器将此噪声滤除，从而恢复原信号。

(D) 发现此信号夹杂了一个正弦噪声noise= cos(2π*60t)，请设计一个高通滤波器将此噪声滤除，从而恢复原信号。

要求：
（1）请写出具体的MATLAB程序，并详细解释每条程序（2）画出滤波前后信号的频谱图
（3）画出所设计滤波器的幅频和相频特性图，并写出具体参数
参数计算：
根据题目要求，开始选取Wp=2*60π，Ws=2*140π。

后来经老师指点，为了将阻带里的信号更好的滤除，通带里的信号更好的保持，达到较好的滤波效果，通带截止频率选取：Wp=2*70π>2*60π，阻带截止频率选取：Ws=2*120π
<2*140π，输入信号为：
y=cos(2π*60t)+1.2cos(2π*140t)+2sin(2π*220t)
+1.5sin(2π*300t) 可知信号最高频率为2*300*π/(2
π)=300Hz。

由奈奎斯特抽样定理得，fs>=2*300=600(Hz)，这里为了得到更好的抽样效果，同时简化计算，选取
fs=1000Hz。

下面计算关于π的归一化频率：
通带截止频率：wp=Wp/fs=0.14*π
阻带截止频率：ws=Ws/fs=0.24*π
软件介绍：
简介：MATLAB 是一种用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。

使用 MATLAB，您可以较使用传统的编程语言（如 C、C++ 和 Fortran）更快地解决技术计算问题。

MATLAB 的应用范围非常广，包括信号和图像处理、通讯、控制系统设计、测试和测量、财务建模和分析以及计算生物学等众多应用领域。

附加的工具箱（单独提供的专用 MATLAB 函数集）扩展了MATLAB 环境，以解决这些应用领域内特定类型的问题。

MATLAB 提供了很多用于记录和分享工作成果的功能。

可以将您的 MATLAB 代码与其他语言和应用程序集成，来分发您的 MATLAB 算法和应用。

主要功能:
1.此高级语言可用于技术计算
2.此开发环境可对代码、文件和数据进行管理
3.交互式工具可以按迭代的方式探查、设计及求解问题
4.数学函数可用于线性代数、统计、傅立叶分析、筛选、优化以及数值积分等
5.二维和三维图形函数可用于可视化数据
6.各种工具可用于构建自定义的图形用户界面
7.各种函数可将基于 MATLAB 的算法与外部应用程序和语言（如 C、C++、Fortran、Java、COM 以及 Microsoft Excel）集成
题目分析：
根据题目要求，综合在《数字信号处理》中所学的知识以及老师的建议与讲解，此设计用FIR、IIR都可实现。

利用MATLAB中的专用函数进行编写，最终确定设计方案如下：
1.窗函数法设计FIR数字滤波器
2.切比雪夫1型高通滤波器
第一方案：窗函数法设计FIR 数字滤波器 窗函数设计法的基本思想是用FIRDF 逼近希望的滤波器。

设希望逼近的滤波器的频率为()
ωj d e H ，其单位脉冲响应用()n h d 表示，为了设计简单方便，通常选择()
ωj d e H 具有片段常数特性的理想滤波器因此()n h d 是无限长的非因果序列，不能直接作为FIRDF 的单位脉冲响应。

窗函数设计法就是截取()n h d 有限长的一段因果序列，并用适合的窗函数进行加权作为FIRDF 的的单位脉冲响应()n h 。

实验程序：
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 滤波器部分 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% wp=0.24*pi; %关于π的归一化通带截止频率
ws=0.14*pi; %关于π的归一化阻带截止频率
DB=wp-ws; %过渡带宽
N0=ceil(6.2*pi/DB); %计算所需h(n)长度N0，ceil(x)取大于等于x 的最小整数
N=N0+mod(N0+1,2); %确保h(n)长度N 是奇数
wc=(wp+ws)/2/pi; %计算理想高通滤波器通带截止频率（关于π归一化） hn=fir1(N-1,wc,'high',hanning(N));%调用fir1计算高通FIRDF 的h(n) Fs=1000; %抽样频率
[H,w1]=freqz(hn,1,N,Fs);%求滤波器幅度响应，设置最大幅度为1
plot(w1,abs(H)); %画图，滤波器幅度响应
title('滤波器幅度响应'); %设置图像窗口标题
figure(2); %创建图像窗口（2）
freqz(hn); %画图，滤波器幅度响应(db)和相位响应
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 输入信号部分 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% t=0:0.001:1.999; %设置t 变量范围，和步长
n=2000; %抽样点数
Fs=1000; %抽样频率
y=cos(2*pi*60*t)+1.2*cos(2*pi*140*t)+2*sin(2*pi*220*t)+1.5*sin(2*pi*3 00*t);%输入信号
y1=fft(y); %输入信号的傅里叶变换
y2=fftshift(y1); %输入信号傅里叶变换重新排布，使数据与频率对应f=(0:1999)*Fs/n-Fs/2; %计算频率f
hold on; %保持图形
figure(3); %创建图像窗口（3）
plot(f,abs(y2),'b'); %画图，输入信号频谱图
title('输入信号频谱图'); %设置图像窗口标题
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 输出信号部
分 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
G=fftfilt(hn,y); %输入信号y通过滤波器
G1=fft(G); %滤波后，输出信号傅里叶变换
G2=fftshift(G1); %输出信号傅里叶变换重新排布，使数据与频率对应figure(4); %创建图像窗口（4）
plot(f,abs(G2));grid %画图，输出信号频谱图
title('输出信号频谱图'); %设置图像窗口标题
实验图形：
滤波器的幅频特性，可以看出，再阻带截止频率70Hz以下的幅频响应基本为0，通带截止频率120Hz以上的幅频响应基本为1。

系统滤波器有较好的滤波效果。

上图幅频响应(db)通带较好保持，阻带有较大衰减。

下图相频响应，在通带内滤波器有良好的线性相位。

输入信号的频谱图（含有低频部分60Hz）
输出信号，可以看出低频部分60Hz已被滤除
第二方案：雪比切夫1型
cheblap,cheblord,cheby1是切比雪夫I 型滤波器设计函数。

其调用格式如下： ①．［z,p,G ］=cheblap(N,Rp)
计算N 阶归一化系统函数的零、极点和增益因子。

返回长度为N 的列向量z 和p 分别给出N 个零点和N 个极点的位置，G 表示滤波器增益。

得到的系统函数如下：
))
())...(2())(1(())())...(2())(1(()()()(N p s p s p s N z s z s z s G s D s P s H a a a ------== 式中，z(k)和p(k)分别为向量z 和p 的第k 个元素。

②．[N,wpo]=cheblord(wp,ws,Rp,As)
计算滤波器的阶数和通带截止频率，wp 和ws 为数字滤波器的通带边界频率和阻带边界频率的归一化值，Rp 和As 为通带最大衰减和阻带最小衰减。

③．[B,A]=cheby1(N,Rp,wpo,'ftype')
N,wpo 分别为滤波器阶数和通带截止频率，该式计算系统函数分子和分母系数多项式的系数向量B 和A 。

N N N
N z
N A z N A z A A z N B z N B z B B z A z B z H --------++++++++++==)1()(...)2()1()1()(...)2()1()()()()1(1)1(1 式中，B(k)和A(k)分别为向量B 和A 的第k 个元素。

自定义函数（freqz_m ）：
function [db,mag,pha,grd,w] = freqz_m(b,a);
% Modified version of freqz subroutine
% ------------------------------------
% [db,mag,pha,grd,w] = freqz_m(b,a);
% db = Relative magnitude in dB computed over 0 to pi radians % mag = absolute magnitude computed over 0 to pi radians % pha = Phase response in radians over 0 to pi radians % grd = Group delay over 0 to pi radians
% w = 501 frequency samples between 0 to pi radians % b = numerator polynomial of H(z) (for FIR: b=h) % a = denominator polynomial of H(z) (for FIR: a=[1]) %
[H,w] = freqz(b,a,1000,'whole');
H = (H(1:1:501))'; w = (w(1:1:501))';
mag = abs(H);
db = 20*log10((mag+eps)/max(mag));
pha = angle(H);
% pha = unwrap(angle(H));
grd = grpdelay(b,a,w);
% grd = diff(pha);
% grd = [grd(1) grd];
% grd = [0 grd(1:1:500); grd; grd(2:1:501) 0];
% grd = median(grd)*500/pi;
实验程序：%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 滤波器部分 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% wp1=0.24*pi; %关于π的归一化通带截止频率
ws1=0.14*pi; %关于π的归一化阻带截止频率
Rp=0.1; %通带最大衰减
Ar=15; %阻带最小衰减
[N1,wn1]=cheb1ord(wp1/pi,ws1/pi,Rp,Ar);%在给定滤波器性能的情况下，选择
雪比切夫1型滤波器的阶数和截止频率Wn1
[b1,a1]=cheby1(N1,Rp,wn1,'high'); %根据阶数，通带最大衰减，截止频率，设置雪比切夫1型高通滤波器
[db1,mag1,pha1,grd1,w1]=freqz_m(b1,a1);%求滤波器幅度响应（db），幅度响应，相位响应，群延迟
figure(1) %创建图像窗口（1）
subplot(3,4,1); %分割图像窗口三行四列，第一小窗口
plot(w1/pi,mag1); %画图，幅度响应
axis([0 1 0 1.1]); %axis([xmin xmax ymin ymax])设置x轴和y轴的表
示范围
ylabel('(a)'); %设置纵坐标y的名称
title('幅度响应'); %设置小窗口标题
subplot(3,4,2); %分割图像窗口三行四列，第二小窗口
plot(w1/pi,db1); %画图，幅度响应（db）=20log（幅度响应）
axis([0 1 -40 5]); %axis([xmin xmax ymin ymax])设置x轴和y轴的表
示范围
title('幅度响应（db）'); %设置小窗口标题
subplot(3,4,3); %分割图像窗口三行四列，第三小窗口
plot(w1/pi,pha1/pi); %画图，相位响应
axis([0 1 -1 1]); %axis([xmin xmax ymin ymax])设置x轴和y轴的表
示范围
title('相位响应'); %设置小窗口标题
subplot(3,4,4); %分割图像窗口三行四列，第四小窗口
plot(w1/pi,grd1); %画图，群延迟
axis([0 1 0 10]); %axis([xmin xmax ymin ymax])设置x轴和y轴的表
示范围
title('群延迟'); %设置小窗口标题
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 输入信号部分 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t=0:0.001:1.999; %置t变量范围，和步长
n=2000; %抽样点数
Fs=1000; %抽样频率
y=cos(2*pi*60*t)+1.2*cos(2*pi*140*t)+2*sin(2*pi*220*t)+1.5*sin(2*pi*3
00*t);%输入信号
y1=fft(y); %输入信号的傅里叶变换
y2=fftshift(y1); %输入信号傅里叶变换重新排布，使数据与频率对应
f=(0:1999)*Fs/n-Fs/2; %计算频率f
hold on; %保持图形
figure(2); %创建图像窗口（2）
plot(f,abs(y2),'b');grid%画图，输入信号频谱图
title('输入信号频谱图'); %设置图像窗口标题
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 输出信号部
分 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
G=filter(b1,a1,y); %输入信号y通过滤波器
G1=fft(G); %滤波后，输出信号傅里叶变换
G2=fftshift(G1); %输出信号傅里叶变换重新排布，使数据与频率对应
figure(3); %创建图像窗口（3）
plot(f,abs(G2));grid %画图，输出信号频谱图
title('输出信号频谱图'); %设置图像窗口标
滤波器特性：
幅频响应及幅频响应(db)，通带较好保持，阻带有较大衰减。

输入信号的频谱图（含有低频部分60Hz）滤波后频谱图：
输出信号，可以看出低频部分60Hz已被滤除
实验心得：
由于，数字信号处理是上个学期修的课程，所以开始的时候，内容有点忘了。

老师布置课程设计任务以后，把书本重新复习了一遍，对内容有所回忆。

也大致了解了课程设计的内容和任务。

然后按照数字信号处理知识，将课程设计的任务参数算出来。

根据参数选择设计方案。

之后，用matlab加以实现。

这里遇到另一个问题，我们对matlab这个软件不太熟悉。

语法也不熟悉，我们通过对照书本学习，又到图书馆借了一些matlab关于数字信号处理的资料，上网查一些相关函数
和资料。

最后一步一步将滤波器，信号的频谱图，信号如何通过滤波器（滤波器是怎样对信号进行滤波的），都实现了。

设计过程中，我学到了一个思想，在确认原理基本正确的时候，如果程序还存在问题，有可能是语法问题，这时候应该上网查找资料，或者查找其他书籍资料，以快速、有效解决。

比如FIR滤波器，我们在做信号如何通过滤波器时，试图把h(n)快速傅里叶变（fft）换后H，与信号y fft后的Y直接相乘，然后输出信号G=H*Y，这是课本原理，在matlab上语法不是这样的，他有专门滤波函数G=fftfilt(hn,y)。

通过这次课程设计，我们对《数字信号处理》这门课程有了更深一些的了解。

对matlab也有了一定的认识和了解，熟悉了matlab在数字信号处理中的一些基本语法。

参考资料：
[1] 《数字信号处理——原理、实现及应用》.高西全，丁玉美，阔永红.北京：电子工业出版社，2008.
[2] 《数字信号处理实践教程》.杨述斌，李永全.武汉：华中科技大学出版社，2007.
[3] 《数字信号处理及MATLAB实现（第二版）》.余成波，陶红艳，杨菁，杨如民.北京：清华大学出版社，2008.
[4]《MATLAB在数字信号处理中的应用（第2版）》.薛年喜.北京：清华大学出版社，2008.
[5]《数字信号处理基础及MATLAB实现》.周辉，董正宏.北京：北京希望电子出版社，中国林业出版社.2006.
光信息08-3班张静2011年7月1日。