MATLAB FIR 滤波器.doc

Matlab课程设计报告目录摘要21.数字滤波器与MATLAB概述41.1数字滤波器41.2 MATLAB简介52.窗函数72.1用窗函数设计FIR数字滤波器的基本方法72.2典型窗函数及其调用格式82.2.1典型窗函数82.2.2典型窗函数的调用格式103.基于窗函数的FIR滤波器的MATLAB实现113.1 理想低通滤波器实现113.2系统各响应函数113.3 滤波器主函数123.4 滤波器主程序解析154.实验结果及其分析174.1原信号的时域波形和处理后时域波形174.2理想带通与实际带通的脉冲响应174.3各窗函数幅频相频波形184.4所构造带通滤波器幅频相频波形204.5 原信号频域波形与处理后信号的频域波形235.收获体会246.参考文献25摘要随着信息时代和数字世界的到来，数字信号处理已成为今一门极其重要的学科和技术领域。数字信号处理在通信、语音、图像、自动控制、雷达、军事、航空航天、医疗和家用电器等众多领域得到了广泛的应用。在数字信号处理应用中，数字滤波器十分重要并已获得广泛应用。数字滤波是数字信号处理的重要内容，数字滤波器可分为IIR和FIR两大类。对于FIR数字滤波器的设计，可以根据所给定的频率特性直接设计，文中采用的设计方法是窗函数法。本文根据FIR滤波器的特点，在MATLAB坏境下用窗函数设计FIR数字滤波器，并对信号进行分析，最后给出了FIR带通滤波器对信号滤波的效果。关键词 数字滤波器；FIR ；MATLAB ；窗函数Abstract With the information age and the advent of the digital world, digital signal processing has become today a very important disciplines and technical fields. Digital signal processing in communications, voice, image, automatic control, radar, military, aerospace, medical and household appliances and many other fields has been widely used. In digital signal processing applications, the digital filter is very important and has been widely used. Digital filtering is an important part of digital signal processing, digital IIR and FIR filters can be divided into two broad categories. For the FIR digital filter design, can be given according to the frequency characteristics of the direct design, design method used in the text its a window function method. Based on the characteristics of FIR filters in MATLAB under the bad with the window function throughout the design FIR digital filters, and signal analysis, and finally gives a FIR band-pass filter the signal filtering effect. Key words ： digital filter; FIR; MATLAB; window function 1.数字滤波器与MATLAB概述 1.1数字滤波器数字滤波器是指完成信号滤波处理功能的，用有限精度算法实现的离散时间线性非时变系统，其输入是一组数字量，其输出是经过变换的另一组数字量。因此，数字滤波器本身既可以是用数字硬件装配成的一台完成给定运算的专用的数字计算机，也可以将所需要的运算编成程序，让通用计算机来执行。数字滤波器，输入输出均为数字信号，通过一定的运算关系，改变输入信号中所含频率成分的相对比例，或则滤除某些频率成分的器件。对于数字滤波器而言，若系统函数为H(z),其脉冲响应为h(n),输入时间序列为x(n),则它们在时域内的关系式如下：y(n)=h(n)*x(n) (公式1-1)在Z域内，输入和输出存在如下关系： Y(z)= H(z)X(z) (公式1-2)式中，X(z)、Y(z)分别为x(n)和y(n)的Z变换。在频域内，输入和输出则存在如下关系：Y(j)=H(j)X(j) (公式1-3)式中，H(j)是数字滤波器的频率特性；X(j)Y(j)分别为x(n)和y(n)的频谱，而为数字角频率。数字滤波器可以有很多种分类方法，但总体上可分为两大类。一类称为经典滤波器，即一般的滤波器，其特点是输入信号中的有用成分和希望滤除的成分占用不同的频带，通过合适的选频滤波器可以实现滤波。例如，若输入信号中有干扰，信号和干扰的频带互不重叠，则可滤出信号中的干扰得到纯信号。但是，如果输入信号中信号和干扰的频带相重叠，则干扰就不能被有效的滤出。另一类称为现代滤波器，如维纳滤波器、卡尔曼滤波器等，其输入信号中有用信号和希望滤除的成分频带重叠。对于经典滤波器，从频域上也可以分为低通、高通、带通和带阻滤波器。从时域特性上看，数字滤波器还可以分为有限冲激响应数字滤波器（FIR）和无限冲激响应数字滤波器（IIR）。对于有限冲激响应数字滤波器（FIR），其输出y(n)只取决于有限个过去和现在的输入，x(n),x(n-1),x(n-m),滤波器的输入输出关系可表示为 y(n)= (公式1-4)对于无限冲激响应数字滤波器（IIR），它的输出不仅取决于过去和现在的输入，而且还取决于过去的输出，其差分方程为y(n)+= (公式1-5)该差分方程的单位冲激响应是无限延续的。1.2 MATLAB简介MATLAB是美国MathWorks公司开发的一种功能极其强大的高技术计算语言和内容极其丰富的软件库，集数值计算、矩阵运算和信号处理与显示于一身。该软件最初是由美国教授Cleve Moler创立的。1980年前后，他在教线性代数课程时，发现用其他高级语言编程时极不方便，便构思开发了MATLAB，即矩阵实验室（Matrix Laboratory）。该软件利用了当时代表数值线性代数领域最高水平的EISPACK和LINPACK两大软件包，并且利用Fortran语言编写了最初的一套交互式软件系统，MATLAB的最初版本便由此产生了。 最初的MATLAB由于语言单一，只能进行矩阵的运算，绘图也只能用原始的描点法，内部函数只有几十个，因此功能十分简单。1984年该公司推出了第一个MATLAB的商业版，并用C语言作出了全部改写。现在的MATLAB程序是MathWorks公司用C语言开发的，第一版由steve Bangert主持开发编译解释程序，Steve Kleiman完成图形功能的设计，John Little和Cleve Moler主持开发了各类数学分分析的子模块，撰写用户指南和大部分的M文件。接着又添加了丰富的图形图像处理、多媒体功能、符号运算和与其它流行软件的接口功能，使MATLAB的功能越来越强大。MTALAB系统主要由以下五个部分组成：(1)MATALB语言体系。 MATLAB是高层次的矩阵数组语言，具有条件控制、函数调用、数据结构、输入输出、面向对象等程序语言特性。利用它既可以进行小规模端程，完成算法设计和算法实验的基本任务，也可以进行大规模编程，开发复杂的应用程序。(2)MATLAB工作环境 。这是对MATLAB提供给用户使用的管理功能的总称。包括管理工作空间中的变量据输入输出的方式和方法，以及开发、调试、管理M文件的各种工具。(3)图形句相系统 。这是MATLAB图形系统的基础，包括完成2D和3D数据图示、图像处理、动画生成、图形显示等功能的高层MATLAB命令，也包括用户对图形图像等对象进行特性控制的低层MATLAB命令，以及开发GUI应用程序的各种工具。(4)MATLAB数学函数库。这是对MATLAB使用的各种数学算法的总称。包括各种初等函数的算法，也包括矩阵运算、矩阵分析等高层次数学算法。(5)MATLAB应用程序接口(API)。这是MATLAB为用户提供的一个函数库，使得用户能够在MATLAB环境中使用C程序或FORTRAN程序，包括从MATLAB中调用于程序(动态链接)，读写MAT文件的功能。 除此之外，MATLAB系统还具有如下特点：（1）具有易学易用的语言体系；（2）具有交互式的工作环境；（3）具有多层面的图像处理系统；（4）具有丰富高效的MATLAB工具箱；（5）具有便利的程序接口（API）;（6）应用领域广泛； （7）嵌入了面向对象编程语言。2.窗函数2.1用窗函数设计FIR数字滤波器的基本方法设计思想：从时域从发，设计逼近理想。设理想滤波器的单位脉冲响应为。以低通线性相位FIR数字滤波器为例。 (公式2-1)一般是无限长的，且是非因果的，不能直接作为FIR滤波器的单位脉冲响应。要想得到一个因果的有限长的滤波器h(n)，最直接的方法是截断，即截取为有限长因果序列，并用合适的窗函数进行加权作为FIR滤波器的单位脉冲响应。按照线性相位滤波器的要求，h(n)必须是偶对称的。对称中心必须等于滤波器的延时常数，即 (公式2-2)用矩形窗设计的FIR低通滤波器，所设计滤波器的幅度函数在通带和阻带都呈现出振荡现象，且最大波纹大约为幅度的9%，这个现象称为吉布斯效应。根据过渡带宽及阻带衰减要求，选择窗函数的类型并估计窗口长度N（或阶数M=N-1），窗函数类型可根据最小阻带衰减As独立选择，因为窗口长度N对最小阻带衰减As没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度N，设待求滤波器的过渡带宽为w，它与窗口长度N近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的N，在N和窗函数类型确定后，即可调用MATLAB中的窗函数求出窗函数wd（n）。根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n)，如果给出待求滤波器频率应为Hd，则理想的单位脉冲响应可以用下面的傅里叶反变换式求出： (公式2-3) 在一般情况下，hd(n)是不能用封闭公式表示的，需要采用数值方法表示；从w=0到w=2采样N点，采用离散傅里叶反变换（IDFT）即可求出。用窗函数wd(n)将hd(n)截断，并进行加权处理，得到 (公式2-4)如果要求线性相位特性， 则h(n)还必须满足： (公式2-5)根据上式中的正、 负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类。 例如， 要设计线性相位低通特性可选择h(n)=h(N-1-n)一类，而不能选h(n)=-h(N-1-n)一类。 验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用freqz子程序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。2.2典型窗函数及其调用格式2.2.1典型窗函数1、矩形窗(Rectangle Window) (公式2-6)其频率响应和幅度响应分别为：， (公式2-7)2、三角形窗(Bartlett Window) (公式2-8)其频率响应为： (公式2-9)3、汉宁(Hanning)窗，又称升余弦窗 (公式2-10)其频率响应和幅度响应分别为： (公式2-11)4、汉明(Hamming)窗，又称改进的升余弦窗 (公式2-12)其幅度响应为： (公式2-13) 5、布莱克曼(Blankman)窗，又称二阶升余弦窗 (公式2-14)其幅度响应为： (公式2-15)6、凯泽(Kaiser)窗 (公式2-16)其中：是一个可选参数，用来选择主瓣宽度和旁瓣衰减之间的交换关系，一般说来，越大，过渡带越宽，阻带越小衰减也越大。I0()是第一类修正零阶贝塞尔函数。若阻带最小衰减表示为，的确定可采用下述经验公式： (公式2-17)若滤波器通带和阻带波纹相等即p=s时，滤波器节数可通过下式确定： (公式2-18)其中： (公式2-19)在MATLAB中，实现矩形窗的函数为boxcar和rectwin，其调用格式如下： w=boxcar(N) w=rectwin(N)其中N是窗函数的长度，返回值w是一个N阶的向量，它的元素由窗函数的值组成。实际上，w=boxcar(N)等价于w=ones(N,1)。2.2.2典型窗函数的调用格式1. 在MATLAB中，实现三角窗的函数为boxcar，调用格式为： w=boxcar(N)2. 在MATLAB中，实现三角窗的函数为triang，调用格式为： w=triang(N)3. 在MATLAB中，实现汉宁窗的函数为hann，调用格式如下： w=hann(N)w=hann(N, sflag)Hann函数中的参数sflag为采样方式，其值可取symmetric(默认值)或periodic。当sflagsymmetric时，为对称采样；当sflagperiodic时，为周期采样，此时hann函数计算N+1个点的窗，但是仅返回前N个点。4. 在MATLAB中，实现海明窗的函数为hamming，调用格式分别如下： w=hamming (N)w=hamming (N,sflag)其中sflag的用法同上。5. 在MATLAB中，实现布拉克曼窗的函数为blackman，调用格式如下： w=blackman (N)w=blackman (N,sflag)6. 在MATLAB中，实现切比雪夫窗的函数为chebwin，调用格式为： w=chebwin (N,r)其中r 表示切比雪夫窗函数的傅里叶变换旁瓣幅度比主瓣低rdB(其默认值为100dB)，且旁瓣是等纹波的。7. 在MATLAB中，实现巴特里特窗的函数为bartlett，调用格式为： w=bartlett (N)3.基于窗函数的FIR滤波器的MATLAB实现3.1 理想低通滤波器实现 通过先产生一个理想的低通滤波器变换为一理想的带通滤波器,作为实际滤波器的模拟对象.该理想低通滤波器由ideal_lp.m文件产生. 该程序代码如下所示: ideal_lp.m文件:function hd=ideal_lp(w,N);alpha=(N-1)/2;n=0:(N-1);m=n-alpha+eps;hd=sin(w*m)./(pi*m);3.2系统各响应函数 该函数为求取系统的绝对幅度响应、相对的db值幅度响应、相位响应和群延时响应的函数,该函数由Freqz_m.m文件实现. 该程序代码如下所示:Freqz_m.m文件function db,mag,pha,grd,w=freqz_m(b,a)%求取系统的绝对幅度响应、相对的db值幅度响应、相位响应和群延时响应的函数%db为相对振幅(dB)%mag为绝对振幅%pha为相位响应%grd为群延时%w为频率样本点向量%b为Ha(z)分子多项式系数(对FIR而言，b=h)%a为Hz(z)分母多项式系数(对FIR而言，a=1) H,w=freqz(b,a,1000,whole); %freqz显示数字滤波器频域中的图形H=(H(1:501);w=(w(1:501);mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);grd=grpdelay(b,a,w);3.3 滤波器主函数该带通滤波器由Daitong.m文件实现.该程序的具体代码如下:Daitong.m文件Function daitongN1=1024;t=0:1:N1-1;fs=5000;s=(sin(2*100*pi*t/fs)+sin(2*pi*1000*t/fs)+2*sin(2*pi*2000*t/fs)+randn(1,length(t/fs); %处理前的信号函数Y=fft(s); %原信号的傅里叶变换ws1 = 0.2*pi; wp1 = 0.35*pi;wp2 = 0.65*pi; ws2 = 0.8*pi;As = 60;tr_width = min(wp1-ws1),(ws2-wp2)M = ceil(11*pi/tr_width) + 1 n=0:1:M-1;wc1 = (ws1+wp1)/2; wc2 = (wp2+ws2)/2;hd = ideal_lp(wc2,M) - ideal_lp(wc1,M);%理想带通滤波器w_bla = (blackman(M); %调用窗函数h = hd .* w_bla; %加窗db,mag,pha,grd,w = freqz_m(h,1);%求取系统的绝对幅度响应、相对的db值幅度响应、相位响应和群延时响应的函数delta_w = 2*pi/1000;Rp = -min(db(wp1/delta_w+1:1:wp2/delta_w) As = -round(max(db(ws2/delta_w+1:1:501) Ts=t(2)-t(1); Ws=2*pi/Ts;Wn=Ws/2;y1=filter(h,1,s); %求处理后信号时域波形Y1=fft(y1); %求处理后信号的频谱Ya=abs(Y(1:length(t)/2);Ya1=abs(Y1(1:length(t)/2);f=linspace(0,Wn/(2*pi),length(t)/2);figure(1); %打印理想脉冲响应波形和实际脉冲响应的波形subplot(1,2,1);stem(n,hd); title(理想脉冲响应)axis(0 M-1 -0.4 0.5); ylabel(hd(n)subplot(1,2,2); stem(n,w_bla);title(Blackman Window)axis(0 M-1 0 1.1); ylabel(w(n)subplot(2,2,3); stem(n,h);title(实际脉冲响应)axis(0 M-1 -0.4 0.5); ylabel(h(n)figure(2); %打印滤波器频率响应和相位响应freqz(h,1);grid;ylabel(Decibels)axis(0 1 -150 10);figure(3); %打印原信号时域波形和处理后的时域波形subplot(1,2,1)plot(t,s);title(原信号时域图);axis(0,300,-5,5);grid;ss=ifft(Y1(1:length(t);subplot(122)plot(t,ss);title(滤波后的时域图);axis(0,300,-5,5);grid;figure(4); %打印原信号频谱和处理后信号的频谱subplot(1,2,1);plot(f*fs,Ya);axis(0,3000,0,1000);grid;xlabel(f);ylabel(频谱幅值F（w）);title(原信号的频谱);subplot(1,2,2);plot(f*fs,Ya1);axis(0,3000,0,1000);grid;xlabel(f);ylabel(频谱幅值F（w）);title(经过FIR带通滤波器后的信号频谱);3.4 滤波器主程序解析1.需处理信号的产生:s=(sin(2*100*pi*t/fs)+sin(2*pi*1000*t/fs)+2*sin(2*pi*2000*t/fs)+randn(1,length(t/fs);Y=fft(s);如以上程序段,首先定义函数S,函数S为三个频率的正弦信号和一随机信号的叠加生成的信号,作为要处理的信号.然后对该信号做离散傅里叶变换,得到该信号的频域波形函数Y.2.理想滤波器的生成ws1 = 0.2*pi; wp1 = 0.35*pi;wp2 = 0.65*pi; ws2 = 0.8*pi;As = 60;tr_width = min(wp1-ws1),(ws2-wp2)M = ceil(11*pi/tr_width) + 1 n=0:1:M-1;wc1 = (ws1+wp1)/2; wc2 = (wp2+ws2)/2;hd = ideal_lp(wc2,M) - ideal_lp(wc1,M);如上述程序段,首先通过所给的参数,计算出该理想滤波器的通带上下截止频率,然后通过事先生成的理想低通滤波器构造理想带通滤波器hd.3.窗函数的调用w_bla = (boxcar(M); %矩形窗w_bla = (triang(M); %三角窗w_bla = (bartlett(M); %巴特利窗w_bla = (hanning(M); %汗宁窗w_bla = (haimming(M); %海明窗w_bla = (blackman(M); %布莱克曼窗通过上述语句产生需要调用的窗函数.4.加窗h = hd .* w_bla;通过理想带通滤波器的函数与窗函数相乘,得到函数h,通过h模拟理想带通滤波器.5.滤波后信号y1=filter(h,1,s);Y1=fft(y1);通过函数filter得到信号s通过带通滤波器后的信号的时域波形y1.调用函数fft(),得到处理后信号的频域波形.6.波形的打印figure(1);subplot(2,2,1);stem(n,hd); title(Ideal Impulse Response)axis(0 M-1 -0.4 0.5); ylabel(hd(n) subplot(122);plot(f*fs,Ya1);axis(0,3000,0,1000);grid;xlabel(f);ylabel(频谱幅值F（w）);title(经过FIR带通滤波器后的信号频谱);通过调用plot()函数,将原信号时域频域波形,滤波器幅频特性,相频特性,滤波后信号的时域频域波形等打印.4.实验结果及其分析4.1原信号的时域波形和处理后时域波形图5-1 原信号与滤波后的时域信号 该信号为基于布莱克曼窗的FIR带通滤波器处理后的结果.可见高频分量和低频分量减少.4.2理想带通与实际带通的脉冲响应图5-2 理想带通与实际带通的脉冲响应该系统为基于基于布莱克曼窗的FIR带通滤波器的脉冲响应.4.3各窗函数幅频相频波形1.矩形窗幅频相频波形图5-3 矩形窗幅频相频波形2.三角窗幅频相频波形图5-4 三角窗幅频相频波形3.巴特利窗幅频相频波形图5-5 巴特利窗幅频相频波形4.汉宁窗幅频相频波形图5-6 汉宁窗幅频相频波形5.海明窗幅频相频波形图5-7 海明窗幅频相频波形6.布