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

课 程 设 计20011 年 7月 1日 设计题目 学 号专业班级 指导教师 学生姓名 张腾达 吴晔 陈丽娟 杨蕾通信电子电路课程设计 ——数字滤波器的设计 张静 20080302 光信息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);%求滤波器幅度响应，设置最大幅度为1plot(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; %计算频率fhold 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 NN zN 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*300*t);%输入信号y1=fft(y); %输入信号的傅里叶变换y2=fftshift(y1); %输入信号傅里叶变换重新排布，使数据与频率对应f=(0:1999)*Fs/n-Fs/2; %计算频率fhold 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)，通带较好保持，阻带有较大衰减。