1、- 1 -目 录实验 1 离散时间信号的频域分析-2实验 2 FFT 算法与应用-7实验 3 IIR 数字滤波器的设计-12实验 4 FIR 数字滤波器的设计-17- 2 -实验 1 离散时间信号的频域分析一实验目的信号的频域分析是信号处理中一种有效的工具。在离散信号的时域分析中，通常将信号表示成单位采样序列 （n）的线性组合，而在频域中，将信号表示成复变量 e 或 e 的线性组合。通过这样的表示，可以将时域的离散njNj2序列映射到频域以便于进一步的处理。在本实验中，将学习利用 MATLAB 计算离散时间信号的 DTFT 和 DFT,并加深对其相互关系的理解。2、实验原理(1)DTFT 和
2、DFT 的定义及其相互关系。序列 x(n)DTFT 定义为 = ()jwXe它是关于自变量 的复函数，且是以 2 为周期的连续函数。()nxejn可以表示为 ,其中， 和 分别jwX()()()jwjwjwreimXX()jwre()jwim是 实部和虚部；还可以表示为 = ，其中，()je j ()|jjX和 分别是 的幅度函数和相位函数；它们都|jw()arg()j()jw是 的实函数，也是以 2 为周期的周期函数。 序列 的 N 点 DFT 定义为 ，()xn221100()()NNjkjknknXexxW是周期为 N 的序列。 与 的关系： 是对 ）在一个周()XkjeXje期中的谱的
3、等间隔 N 点采样，即 ，而2k|()()|jwNk可以通过对 内插获得，即 ()jXe()Xk- 3 -1 (2/)(1)/202)sin()()iN jkNjwk kXeXeN (2)使用到的 MATLAB 命令有基于 DTFT 离散时间信号分析函数以及求解序列的 DFT 函数。1）基于 DTFT 离散时间信号分析函数有：freqz,real,imag,abs,angle。函数 freqz 可以用来计算一个以 e 的有理分式形式给出的序列的 DTFT 值。jfreqz 的形式多样，常见的有 H=freqz(num,den,w)，其中 num 表示序列有理分式 DTFT 的分子多项式系数，d
4、en 表示分母多项式系数（均按 z 的降幂排列），矢量 w 表示在 02 中给定的一系列频率点集合。freqz 函数的其他形式参见帮助文件。在求出 DTFT 值后，可以使用函数 real,imag,abs 和 angle 分别求出并绘出起实部、虚部和相位谱。例如 ）= 利用函()jXe )8109.56.1290jjee数 freqz 计算出 ，然后利用函数 abs 和 angle 分别求出幅频特性与相位()jHe特性最后利用 plot 命令绘出曲线。2）求解序列 DFT 的函数有：fft,ifft。函数 fft（x）可以计算 R 点序列的 R 点 DFT 值；而 fft（x，N）则计算 R
5、点序列的 N 点 DFT，若 RN，则直接截取 R 点 DFT 的前 N 点，若 RN，则直接截取R点DFT的前N点，若RN，则x先进行补零扩展为N点序列再求N点DFT。函数ifft（X）可以计算R点的谱序列的R点IDFT值；而ifft（X,N）同fft（x，N）的情况。函数conv（x1，x2）可以计算两序列的线性卷积。三 实验设备计算机，MATLAB软件。四 实验内容编制信号产生子程序及本实验的频谱分析主程序。实验中需要用到的基本信号包括：（1）三角波序列： x(n)=n+13,0 13-n ,4 0,其它3n7n程序如下：n2=0:7;x2=13,14,15,16,9,8,7,6;sub
6、plot(221);stem(n2,x2);title(x2序列);grid on;k1=0:7;y21=fft(x2,8);magy21=abs(y21);subplot(222);stem(k1,magy21);title(x2的8点FFT);grid on;k2=0:15;y22=fft(x2,16);magy22=abs(y22);- 9 -subplot(224);stem(k2,magy22);title(x2的16点FFT);grid on;波形如下：图2-1 三角波序列以及它的8点和16点 FFT（2）利用 DFT 计算下面两序列的线性卷积 x(n)=3,-4,6,0,6,-4
7、,h(n)=1,-4,3,4程序如下：clear all; close all; clc; N=9;g=3 -4 6 0 6 -4; h=1 -4 3 4; x=conv(g,h) gf=fft(g,N); hf=fft(h,N); xx=gf.*hf X=ifft(xx,N); grid on;- 10 -波形如下：图2-2 用 DFT 计算的序列 x(n)与 h(n)线性卷积图（3）已知某序列 x（n）在单位圆上的 N=64 等分样点的 z 变换为 X（Z ）=Xk（k）= ，k=0,1,2,63。试分别用 N 点 IFFT 程序和 NNkje/213.0点 FFT 程序计算 ,绘出 及
8、。()()xnIDFTXk: ()xn:|()|k程序如下：k=0:63;N=64;X=1./(1-0.13*exp(-j*2*pi*k/N);x=ifft(X,64)stem(abs(x);xlabel(n);ylabel(x(n);title(25FF.单位圆上等分样点的Z变换图 )grid on;波形如下：- 11 -图2-3 单位圆等分样点Z变换图- 12 -实验 3 IIR 数字滤波器的设计一 实验目的从理论上讲，任何的线性时不变（LTI）离散时间系统都可以看做一个数字滤波器，因此设计数字滤波器实际就是设计离散时间系统。数字滤波器包括 IIR（无限冲击响应）和 FIR（有限冲击响应）
9、型，在设计时通常采用不同的方法。本实验通过使用 MATLAB 函数对数字滤波器进行设计和实现，要求掌握IIR 数字巴特沃斯滤波器,数字切比雪夫滤波器的设计原理，设计方法和设计步骤；能根据给定的滤波器指标进行滤波器设计；同时也加深学生对数字滤波器的常用指标和设计过程的理解。二 实验原理在 IIR 滤波器的设计中，常用的方法是：先根据设计要求寻找一个合适的模拟原型滤波器 Ha(s),然后根据一定的准则将比模拟原型滤波器转换为数字滤波器 Hd(z),即我们需要设计的数字滤波器。IIR 滤波器的阶数就等于所选的模拟原型滤波器的阶数，所以其阶数确定主要是在模拟原型滤波器设计中进行的。IIR 数字滤波器的
10、设计方法如下：（1）冲击响应不变法。对满足设计要求的模拟原型滤波器 Ha(s)进行部分分式展开为;基于 = ，可得1()NkkAHass()dnh)aT11()NkdsTkHzezA(2)双线性变换法。对设计要求中给出的边界频率进行预畸处理，然后用得到的频率进行模拟滤波器设计，得到模拟原型滤波器 Ha(s);用双线性变换法求出数字滤波器 。1()()azdzHs一般来说，在要求时域冲激响应能模仿模拟滤波器的场合，一般使用冲激响应不变法。冲击响应不变法一个重要特点是频率坐标的变换是线性的，因此如果模拟滤波器的频率响应带限于折叠频率的话，则通过变换后滤波器的频率响应可不失真地反映原响应与频率的关系
11、。与冲激响应不变法比较，双线性变换的主要优点是靠频率的非线性关系得到 s 平面与 z 平面的单值一一对应关系，整个值对应于单位圆一周。所以从模拟传递函数可直接通过代数置换得到数字滤波器的传递函数。MATLAB 提供了一组标准的数字滤波器设计函数，大大简化了滤波器的设计过程。- 13 -（1） butter。功能：Butterworth 模拟/数字滤波器设计如图；格式：b,a=butter(n,wn,ftype,s),b,a=butter(n,wn,ftype)。（2） Cheby1、cheby2。功能：chebyshev 、chebyshev 型模拟/数字滤波器设计如图所示。格式：b,a=ch
12、eby1(n,Rp,wn,ftype,),b,a=cheby2(n,Rs,wn,ftype)IIR 滤波器设计函数还包括：buttord,chebwin,cheb1ord,cheb2ord,ellip,ellipord 等。还可以用N,Wn=buttord(Wp,Ws,Rp,Rs)的 MATLAB 命名估算一个 Butterworth 滤波器的阶数。三实验设备计算机，MATLAB 软件。四实验内容利用 MATLAB 编程方法或利用 MATLAB 中的 fdatool 工具设计不同功能的IIR 数字滤波器。（1）基于 chebyshev 型模拟滤波器原型使用冲激不变转换方法设计数字滤波器，要求参
13、数为通带截止频率 =0.213 ；通带最大衰减 =1dB;阻带截pWpA止频率 =0.413 ；阻带最小衰减 =35dB。sWsA程序如下：clear;wp=0.213*pi;ws=0.413*pi;Ap=1;As=35;T=1; Rip=10(-Ap/20);Atn=10(-As/20);OmegaP=wp*T;OmegaS=ws*T;n,Wn=cheb1ord(OmegaP,OmegaS,1,35,s); b,a=cheby1(n,1,Wn,low,s); freqs(b,a); bz,az=impinvar(b,a,T);H,W=freqz(bz,az,512,T); - 14 -ma=
14、20*log10(abs(H),pha=20*log10(unwrap(angle(H),hi=impz(bz,az);ni=step(bz,az);subplot(2,2,1),plot(W,ma);title(幅频特性(dB);xlabel(w(/pi);ylabel(dB);subplot(2,2,2),plot(W,pha);title(相频特性);xlabel(w(/pi);ylabel(pha(/pi);subplot(2,2,3),plot(hi);title(单位冲击响应);xlabel(n);ylabel(h(n);subplot(2,2,4),plot(ni); title
15、(单位阶跃响应); xlabel(n);ylabel(h(n);bz,az; 波形如下：图3-1 基于chebyshevI 模型的数字滤波器- 15 -(2)基于 Butterworth 型模拟滤波器原型使用双线性变换方法设计数字滤波器的，要求参数为通带截止频率 =0.24 ；通带最大衰减 =1.3dB;阻带截止pWpA频率 =0.13 ；阻带最小衰减 =40dB。sWsA程序如下：clear;wp=0.4*pi; ws=0.13*pi; Ap=1.3; As=40; fs=1; ts=1/fs; wp2=2*fs*tan(wp/2*ts); ws2=2*fs*tan(ws/2*ts); n,
16、wn=buttord(wp2,ws2,Ap,As,s) z,p,k=buttap(n) bap,aap=zp2tf(z,p,k) b,a=lp2lp(bap,aap,wn) bz,az=bilinear(b,a,fs); h,w=freqz(bz,az); subplot(2,1,1);plot(w/pi,abs(h);grid onxlabel(频率 );ylabel(幅度 );subplot(2,1,2);plot(w/pi,20*log10(abs(h);grid onxlabel(频率 );ylabel(幅度 ); 波形如下:- 16 -图 3-2 基于 Butterworth 模型的
17、数字滤波器- 17 -实验 4 FIR 数字滤波器的设计一实验目的掌握用窗函数法设计 FIR 数字滤波器的原理与其设计步骤；熟悉线性相位数字滤波器的特性。学习编写数字滤波器的设计程序的方法，并能进行正确编程；根据给定的滤波器指标，给出设计步骤。二实验原理如果系统的冲击响应 为已知，则系统的输入-输出关系为:()dhn*()dyx对于低通滤波器，只要设计出低通滤波器的冲击响应函数，就可以由式得到系统的输出了。假设所希望的数字滤波器的频率响应为 H (e )，它是频域djw的周期函数，周期为 2 ，那么它与 H (e )相对应的傅里叶系数为:djw1()(jjndhne但是将 h （n）作为滤波器
18、脉冲响应有两个问题：一是它是无限长的，与dFIR 数字滤脉冲响应有限长这一前提不一致；二是它是非因果的，。对此，要采取两项措施：一是将 截短；二是将其往右()0,d()dhn移。由此得到 的实际频域响应 ，与理想频域响应2()hn120()NjwjwdnHee相近，但不完全一致。理论证明上述现象是对 进行简单截短处理()jwdHe ()d的必然结果，一般称为吉布斯现象，为尽可能的减少吉布斯现象，应对 进()dhn行加窗截取，即以 作为 FIR 滤波器的系数。其中 N 为窗口()()dNhnWn的长度。窗口函数的形状和窗口长度 N 决定了窗函数法设计出的 FIR 滤波()NWn器的性能。设计时，
19、要根据阻带的最小衰减和过渡带宽来选择恰当的窗函数类型和窗长度 N。常用窗函数有矩形窗、海明窗、和布莱克曼窗等。窗函数设计 FIR 滤波器步骤如下：- 18 -（1）给定理想频率响应 的幅频特性和相频特性；()jwdHe（2）求理想单位脉冲响应 ，在实际计算中，可对 采样，并对其求hn()jwdHeIDFT 的 ，用 代替 ；()Mhn()()d（3）根据过渡带宽度和阻带最小衰减，确定窗函数类型和窗长度 N；（4）根据 ,求所需设计的 FIR 滤波器单位脉冲响应；()()dNWn（5）求 ，分析其幅频特性，若不满足要求，可适当改变窗函数形式或jwdHe长度 N，重复上述设计过程，以得到满意的结果
20、。三实验设备计算机，MATLAB 软件。四实验内容(1)用海明窗设计一个 13 阶的 FIR 的低通滤波器，wc=0.13 程序如下：clear all;%海明窗N=13;wc=0.13*pi;h=fir1(N-1,wc/pi,bandpass,hamming(N);h1,wc=freqz(h,1); subplot(2,1,1);plot(wc/pi,20*log10(abs(h1);axis(0,1,-80,10);grid;xlabel(归一化频率/pi);ylabel(幅度/dB);title(N=13,海明窗);波形如下：- 19 -图 4-1 用海明窗设计的低通滤波器(2)用布莱克
21、曼窗设计一个 13 阶的 FIR 的低通滤波器，wc=0.13 程序如下：clear all;%布莱克曼窗N=13;wc=0.13*pi;h=fir1(N-1,wc/pi,bandpass,blackman(N);h1,wc=freqz(h,1); subplot(2,1,1);plot(wc/pi,20*log10(abs(h1);axis(0,1,-80,10);grid;xlabel(归一化频率/pi);ylabel(幅度/dB);title(N=13,布莱克曼窗);波形如下：图 4-2 用布莱克曼窗设计的低通滤波器- 20 -(3)用矩形窗设计线性相位低通滤波器： ，-（ ）0.13,()jwjdHe其余为 0。程序如下：clear all;%矩形窗N=13;wc=0.13*pi;h=fir1(N-1,wc/pi,bandpass,boxcar(N);h1,wc=freqz(h,1); subplot(2,1,1);plot(wc/pi,20*log10(abs(h1);axis(0,1,-80,10);grid;xlabel(归一化频率/pi);ylabel(幅度/dB);title(N=13,矩形窗);波形如下：图 4-3 用矩形窗设计的线性相位低通滤波器
