数字信号处理实验指导书(wcx)

1、数字信号处理数字信号处理 实验指导书实验指导书 王创新王创新 文卉文卉 长沙理工大学电气与信息工程学院电子信息工程教研室长沙理工大学电气与信息工程学院电子信息工程教研室 二零零柒年柒月二零零柒年柒月 第 2 页 共 41 页 目目 录录 前 言 .3 实验一：离散时间信号的时域分析 .4 实验二：离散时间系统的时域分析 .7 实验三：离散时间信号的频域分析 .10 实验四：线性时不变离散时间系统的频域分析 .15 实验五： IIR 数字滤波器的设计 .19 实验六： FIR 数字滤波器的设计.25 附录附录 A MATLAB 系统的常用概念系统的常用概念.30 附录附录 B 信号处理工具箱函数

2、信号处理工具箱函数.35 第 3 页 共 41 页 前前 言言 第 4 页 共 41 页 数字信号处理研究数字序列信号的表示方法，并对信号进行运算，以提取包含 在其中的特殊信息。近几十年来，由于在研究及应用两方面均取得了进展，数字信 号处理领域已日趋成熟。由于计算机的大量使用，从而很容易向学生提供实际环境， 以验证所学的概念和算法。 本指导书编程语言是 MATLAB，它广泛应用于高性能数值计算和可视化。本 书假定读者已具备 MATLAB 基础知识。前面的一些实验帮助学生理解信号处理的 重要概念，后面以设计性实验项目为主，有利于加强对原理的理解并且加强对技术 的应用。附录中给出了本书中用到的 M

3、ATALB 函数及简要解释。 实验一：离散时间信号的时域分析实验一：离散时间信号的时域分析 一实验目的 1、理解时域离散信号的概念，时域离散信号表示方法。 第 5 页 共 41 页 2、熟悉掌握 MATLAB 中产生信号和绘制信号的基本命令。 3、熟悉掌握序列的简单运算，如：加法、标量乘法、时间反转、延时、乘法 等。 二实验相关知识准备 1用到的 MATLAB 命令 运算符： ： . + - * / ; % 基本矩阵： i ones pi rand randn zeros 基本函数： cos exp imag real 数据分析： sum 二维图形： axis grid legend plot

4、 subplot stairs stem title xlable ylable clf 工具箱：sawtooth square 三实验内容 1序列的产生 (1) 程序 11：单位抽样序列的产生和绘制 % Program P1_1 % Generation of a Unit Sample Sequence clf; % Generate a vector from -10 to 20 n = -10:20; % Generate the unit sample sequence u = zeros(1,10) 1 zeros(1,20); % Plot the unit sample seq

5、uence stem(n,u); xlabel(Time index n);ylabel(Amplitude); title(Unit Sample Sequence); axis(-10 20 0 1.2); (2) 程序 12：正弦序列的产生和绘制 第 6 页 共 41 页 % Program P1_2 % Generation of a sinusoidal sequence n = 0:40; f = 0.1; phase = 0; A = 1.5; arg = 2*pi*f*n - phase; x = A*cos(arg); clf;% Clear old graph stem(n

6、,x);% Plot the generated sequence axis(0 40 -2 2); grid; title(Sinusoidal Sequence); xlabel(Time index n); ylabel(Amplitude); axis; 2序列的运算（信号的平滑） 数字信号处理应用的一个常见例子是从被加性噪声污染的信号中移除噪声。假 定信号 sn被噪声 dn所污染，得到一个含有噪声的信号 xn=sn+dn。我们需要 对 xn进行运算，产生一个合理的逼近 sn，对时刻 n 的样本求平均，产生输出信 号是一种简单有效的方法。如：三点滑动平均的信号。 程序 13 实现三点滑

7、动平均的信号运算： yn=(x(n-1)+x(n)+x(n+1)/3 程序如下: % Signal Smoothing by Averaging clf; R = 51; d = 0.8*(rand(R,1) - 0.5); % 产生随机噪声 m = 0:R-1; s = 2*m.*(0.9.m); % 产生未被污染的信号 x = s + d; %产生被噪声污染的信号 subplot(2,1,1); plot(m,d,r-,m,s,g-,m,x,b-.); xlabel(Time index n);ylabel(Amplitude); legend(dn ,sn ,xn ); x1 = 0 0

8、 x;x2 = 0 x 0;x3 = x 0 0; y = (x1 + x2 + x3)/3; subplot(2,1,2); plot(m,y(2:R+1),r-,m,s,g-); legend( yn ,sn ); xlabel(Time index n);ylabel(Amplitude); 第 7 页 共 41 页 3实验问题回答 (1) 命令 clf,axis,title,xlable 和 ylable 的作用是什么？ (2) 程序 12 中正弦序列的频率是多少？怎样可以改变它？哪个参数控制该序 列的振幅？该序列周期是多少？ (3) 程序 13 中加性噪声 dn是什么样的形式？语句

9、x=s+d 代表什么？信号 x1,x2,x3 与信号 x 之间的关系是什么？ 四实验报告要求 1按照实验内容要求完成相关实验程序，并得出相关的实验结果（包括图 形结果） 。 2回答实验中提出的问题。 3总结本次实验结果，按照实验报告格式要求，书写实验报告。 五实验设备 PC 机，MATLAB 软件 实验二：离散时间系统的时域分析实验二：离散时间系统的时域分析 一一实验目的 1、掌握时域离散系统的时域分析方法。 第 8 页 共 41 页 2、通过 MATLAB 仿真一些简单的离散时间系统，并研究它们的时域特性。 二实验相关知识准备 1用到的 MATLAB 命令 语言构造与调试： break en

10、d for if input 基本函数： abs num2str 多项式和内插函数： conv 工具箱： filter impz 其中：y=filter(num,den,x)用于计算线性时不变系统的输出，该命令具有多种形 式，具体使用可参考 MATLAB 的 help 说明。num 和 den 可表示离散时间系统差分 方程中输入与输出变量的系统矩阵。 三实验内容 1线性和非线性系统 例 21 设系统为 yn-0.4yn-1+0.75yn-2=2.2403xn+2.4908xn-1+2.2403xn-2 要求用 MATLAB 程序仿真系统，输入三个不同的输入序列 x1(n),x2(n)和 x(n

11、)=a.x1(n)+b.x2(n)，计算并求出相应的输出响应 y1n,y2n和 yn。 % Generate the input sequences clf; n = 0:40; a = 2;b = -3; x1 = cos(2*pi*0.1*n); x2 = cos(2*pi*0.4*n); x = a*x1 + b*x2; num = 2.2403 2.4908 2.2403; den = 1 -0.4 0.75; ic = 0 0; % Set zero initial conditions y1 = filter(num,den,x1,ic); % Compute the output

12、 y1n y2 = filter(num,den,x2,ic); % Compute the output y2n y = filter(num,den,x,ic); % Compute the output yn yt = a*y1 + b*y2; 第 9 页 共 41 页 d = y - yt; % Compute the difference output dn % Plot the outputs and the difference signal subplot(3,1,1) stem(n,y); ylabel(Amplitude); title(Output Due to Weig

13、hted Input: a cdot x_1n + b cdot x_2n); subplot(3,1,2) stem(n,yt); ylabel(Amplitude); title(Weighted Output: a cdot y_1n + b cdot y_2n); subplot(3,1,3) stem(n,d); xlabel(Time index n);ylabel(Amplitude); title(Difference Signal); 2时不变系统和时变系统 实验内容 22：因果系统的时不变特性研究 编写 MATLAB 程序仿真例 2.1 给出的系统，以产生两个不同的输入序列

14、 x(n) （在例 21 中 x(n)=a.x1(n)+b.x2(n)）和 x(n+D)，计算并画出相应的输出序列 y1n, y2n+D和 y1n-y2n+D 参数：n = 0:40; D = 10; a = 3.0; b = -2; 实验内容23：线性时不变系统的冲激响应计算 设系统为 yn-0.4yn-1+0.75yn-2=2.2403xn+2.4908xn-1+2.2403xn-2 用 y=impz(num,den,N)计算系统的冲激响应的前 N 个样本。 3. 实验问题回答 (1) 修改程序 21 的相关语句，对由加权输入得到的 yn与相同权系数下输出 y1n和 y2n相加得到的 yt

15、n进行比较，这两个序列是否相等？该系统是线 性系统吗？ (2) 比较实验内容 22 中输出序列 yn和 ydn-10，该系统是时不变系统吗？ (3) 计算并画出实验内容 23 的前 40 个样本。 四实验报告要求 1按照实验内容要求完成相关实验程序，并得出相关的实验结果（包括图形 结果） 。 2回答实验中提出的问题。 第 10 页 共 41 页 3. 总结本次实验结果，按照实验报告格式要求，书写实验报告。 五实验设备 PC 机，MATLAB 软件 第 11 页 共 41 页 实验三：离散时间信号的频域分析实验三：离散时间信号的频域分析 一实验目的 1在学习了离散时间信号的时域分析的基础上，对这

16、些信号在频域上进行分 析，从而进一步研究它们的性质。 2熟悉离散时间序列的 3 种表示变换域分析方法：离散时间傅立叶变换 （DTFT）,离散傅立叶变换（DFT）和 Z 变换。 二实验相关知识准备 1用到的 MATLAB 命令 运算符和特殊字符： .* . 语言构造与调试： error function pause 基本函数： angle conj rem 数据分析和傅立叶变换函数： fft ifft max min 工具箱： freqz impz residuez zplane 其中：h=freqz(num,den,w)，可计算离散时间信号(系统)的傅里叶变换，w 为离 散频率点值，一般频率点

17、个数（长度）选为 2 的幂，如 256，512 等。 二三实验内容 1. 离散时间傅里叶变换（DTFT） 定义：()( ) jj n n X ex n e 例 3-1 计算离散时间信号 例 2-1 产生的序列 x(n)的 DTFT 并画出幅度谱和相 位谱。 % Generate discrete time sequences clf; n = 0:40; 第 12 页 共 41 页 a = 2;b = -3; x1 = cos(2*pi*0.1*n); x2 = cos(2*pi*0.4*n); x = a*x1 + b*x2; %计算序列x(n)的DTFT,取频率区间-22 ，点数为128.

18、 w=-2*pi:4*pi/127:2*pi; h=freqz(x,1,w); subplot(2,1,1); plot(w/pi,abs(h); grid; title(|H(ejomega)| 幅度谱); xlabel(omega/pi); subplot(2,1,2) plot(w/pi,angle(h); grid; title(相位谱 argH(ejomega); xlabel(omega/pi); ylabel(以弧度为单位的相位); 2. 离散傅立叶变换 在 MATLAB 中，使用 fft 可以很容易地计算有限长序列 xn的离散傅立叶变 换。此函数有两种形式： y=fft(x)

19、y=fft(x,n) 求出时域信号 x 的离散傅立叶变换，n 为规定的点数，n 的默认值为 所给 x 的长度。当 n 取 2 的整数幂时变换的速度最快。通常取大于又最靠近 x 的幂 次。 （即一般在使用 fft 函数前用 n=2nextpow2(length(x)得到最合适的 n)。 当 x 的长度小于 n 时，fft 函数在 x 的尾部补 0，以构成长为 n 点数据。 当 x 的长度大于 n 时，fft 函数将序列 x 截断，取前 n 点。 一般情况下，fft 求出的函数多为复数，可用 abs 及 angle 分别求其幅度和相位。 观察注意：栅栏效应，截断效应（频谱泄露和谱间干扰） ，混叠失

20、真 例 312： fft 函数最通常的应用是计算信号的频谱。考虑一个由 100hz 和 200hz 正弦信号构成的信号，受零均值随机信号的干扰，数据采样频率为 1000hz。 通过 fft 函数来分析其信号频率成分。 t=0:0.001:1;%采样周期为0.001s，即采样频率为1000hz x=sin(2*pi*100*t)+sin(2*pi*200*t)+1.5*rand(1,length(t);%产生受噪声污染的 正弦波信号 subplot(2,1,1); plot(x(1:50);%画出时域内的信号 第 13 页 共 41 页 y=fft(x,512);%对x进行512点的fft f=

21、1000*(0:256)/512;%设置频率轴（横轴）坐标，1000为采样频率 subplot(2,1,2); plot(f,abs(y(1:257);%画出频域内的信号 实验内容 323：混叠现象与频率分辨率（正确理解奈奎斯特采样定理） 假设现有含有三种频率成分的信号 x(n)=sin(2nf1/fs)+sin(2nf2/fs) +sin(2nf3/fs),其中 f1=20hz,f2=20.5hz,f3=40hz,采样频率 fs=100hz. （1）求 x(n)在 0128 之间的 DFT 的 x(k)； （2）求把（1）中的 x(n)以补 0 方式使其加长到 0512 之间后的 DFT 的

22、 x(k)； （3）求 x(n)在 0512 之间的 DFT 的 x(k)。最后分析三种形式的 DFT 结果图给出 结论。 提示：根据采样定理，在信号采样的过程中，如果采样频率小于信号最高频率 的两倍，那么信号的各次调制频谱就会相互重叠起来，有些频率部分的幅值就与原 始情况不同，因而不能分开和恢复这些部分，取样造成了信号的损失，这种频谱重 叠的出现，就成为“混叠现象” 。 令为采样频率，为信号的最高频率，为了避免混叠现象，则必须使 s f h f ，而频谱分辨率，即在频率轴上所能得到的最小频率间隔 F 又为 hs ff2 N f F s 。所以，在和两个参数中，保持其中一个不变而增加另一个的唯

23、一方法，就 s f h f 是增加在一记录长度内的点数 N。如果与都已给定，则 N 必须满足 s f h f F f N h 2 3. 一维逆快速傅立叶变换 y=ifft(x) y=ifft(x,n) 实验内容 34：频域采样定理的验证 （1） 产生一个三角波序列 x(n) MnMnM Mnn nx 2/ 2/0 )( （2）对 M=40，计算 x(n)的 64 点 DFT，并图示 x(n)和 X(k)=DFTx(n), k=0,1,63。 （3）对（2）中所得 X(k)在0,2pi上进行 32 点抽样得 第 14 页 共 41 页 X1(k)=X(2k), k=0,1,31。并图示。 （4）

24、求 X1（k）的 32 点 IDFT,并图示。即 x1(n)=IDFTX1(k), k=0,1,31 （5）采用周期延拓的方法绘出的波形，评述与 x(n)的关系， 321 )(nx 321 )(nx 并根据频域采样理论加以解释。 4 线线性调频 Z 变换（选做） 离散傅立叶变换（DFT）可以看作信号在 Z 域上沿单位圆的均匀采样。但在实 际应用中，并非整个单位圆上的频谱都有意义。一些情况下，如对于窄带信号，只 希望分析信号所在的一段频带等，采样点的轨迹是一条弧线或圆周。这种需求，就 导致了线性调频 Z 变换（Chirp z 变换）的出现。 Chirp z 变换与 DFT 计算整个频谱的算法不同

25、，它是一种更为灵活的计算频谱 的算法，可以用来计算单位圆上任一段曲线的 Z 变换，作频谱分析时输入的点数 和输出的点数可以不相等，从而达到频域“细化”的目的。 y=czt(x,m,w,a) y=czt(x) 例 35：利用 Chirp z 变换计算滤波器 h（且 h=fir1(30,125/500,boxcar(31)）在 100Hz200Hz 的频率特性，并用图形文件比较 CZT 函数和 FFT 函数。 h=fir1(30,125/500,boxcar(31); fs=1000; f1=100;f2=200; m=1024; y=fft(h,1024); fy=fs*(0:1023)/102

26、4; subplot(2,1,1); plot(fy,abs(y); axis(0,500,0,1.5); w=exp(-j*2*pi*(f2-f1)/(m*fs); a=exp(j*2*pi*f1/fs); z=czt(h,m,w,a);%产生 1024 个 z 值 fz=(f2-f1)*(0:1023)/1024+f1;subplot(2,1,2); plot(fz,abs(z) 5实验问题回答 （1）完成例 31，32 的程序并观察信号序列 x(tn)在时域和频域上分析的 第 15 页 共 41 页 特点。 （2）在实验内容 33 和 34 中按要求完成程序，并回答 33 和 34 中提

27、 出的问题。 (3) 完成例 354 的程序并观察 Chirp z 变换与 fft 变换的特点。 四实验报告要求 1按照实验内容要求完成相关实验程序，并得出相关的实验结果（包括图形 结果） 。 2回答实验中提出的问题。 3总结本次实验结果，按照实验报告格式要求，书写实验报告。 五实验设备 PC 机，MATLAB 软件 第 16 页 共 41 页 实验四：线性时不变离散时间系统的频域分析实验四：线性时不变离散时间系统的频域分析 一实验目的 1掌握离散时间系统的频域表示法。 2在学习了离散时间系统的时域分析的基础上，对这些系统在频域上进行分 析，掌握 LTI 系统的频域分析法。 3掌握熟悉系统函数

28、、传输函数和系统频率响应的计算机仿真方法。? 二实验相关知识准备 1系统传输函数、系统函数 一般称为系统的传输函数，它表征系统的频()( ) jj n n H eh n e () j H e 率响应。|H(ej)|称为幅度响应，而相位H(ej)称为相位响应。 对 h(n)进行 Z 变换，得到 H(Z)，一般称 H(z)为系统的系统响应，系统函数 的一般表示式为： (4-1) 0 0 ( ) ( ) ( ) M i i i N i i i b z Y z H z X z a z 若 H(z)的收敛域包含单位圆，则： (4-2)()( )| j j z e H eH z 12用到的 MATLAB

29、命令 MATLAB 的工具箱，为同一离散系统的不同系统模型间的多种表达方式的变 换，提供了很多功能丰富快捷的 MATLAB 函数，如下表所示： 系统原模型系统转换模型函数 传递函数状态空间Tf2ss 传递函数零极点增益Tf2zp 零极点增益状态空间Zp2ss 零极点增益传递函数Zp2tf 零极点增益二次分式Zp2sos 状态空间传递函数Ss2tf 第 17 页 共 41 页 状态空间零极点增益Ss2zp 状态空间二次分式Ss2sos 二次分式零极点增益Sos2tf 二次分式状态空间Sos2zp 多项式Lattice 结构Sos2ss Lattice 结构多项式Poly2rc 基本函数： abs

30、(计算复数模) angle（计算复数相角） imag(虚部) log10 real(实部) zplane 可以查看零极点 poly2rc 可以直接检测稳定性 三三实验内容 1. 例 41 求由 h(n)=(0.9)nU(n)所表征的系统频率响应 H(ej). 0 + n n=0 ()( )( ) 1 =(0.9) 1 0.91 jj nj n n j j jj H eh n eh n e e e ee clf; w=0:1:500*pi/500; %0 pi axis divided into 501 points H=exp(j*w)./(exp(j*w)-0.9*ones(1,501);

31、magH=abs(H); angH=angle(H); subplot(2,1,1);plot(w/pi,magH);grid; title(Magnitude Response); xlabel(frequency in pi units); ylabel(|H(ej)|); subplot(2,1,2);plot(w/pi,angH);grid; xlabel(frequency in pi units); ylabel(phase in pi radians); title(Phase Response) 2. 例 42 从差分方程求系统频率响应 一个三阶低通滤波器由下面差分方程描述：

32、第 18 页 共 41 页 y(n)=0.0181x(n)-0.0543x(n-1)+0.0543x(n-2)+0.0181x(n-3)+ 1.76y(n-1)-1.1829y(n-2)+0.2781y(n-3) %利用给定的差分方程可求得输入系数 b 和输出系数 a 向量。 clf; b=0.0181 -0.0543 0.0543 0，0181; %coefficient array b a=1 -1.76 1.1829 -0.2781; %coefficient array a m=0:length(b)-1; l=0:length(a)-1; K=500; k=0:1:K; %index

33、 array k for frequencies w=pi*k/K; %0,pi axis divided into 501 points num=b*exp(-j*m*w); %Numerator caculations den = a*exp(-j*l*w); % Denominitor caculations H=num./den; %Frequency response magH= abs(H); angH=angle(H); subplot(2,1,1); plot(w/pi,amgH); grid; axis(0,1,0,1); title(Magnitude Response);

34、 xlabel(frequency in pi units); ylabel(|H(ej)|); subplot(2,1,2);plot(w/pi,angH);grid; xlabel(frequency in pi units); ylabel(phase in pi radians); title(Phase Response);例 41 求截短的理想低通滤波器的冲激响应 理想低通滤波器的冲激响应是无限的，因此实际上不可能实现。在实际应用中， 将冲激响应截短到一个有限值，可以近似地实现该滤波器。然而，将冲激响应截短 后滤波器表现为一个非因果的滤波器。将经过截短后的滤波器的冲激响应向右移动

35、N/2 个样本，可得因果上的近似，如下所示： Nn N n N n nh c 0 , 2 2 sin )( 滤波器的长度为 N+1。 下面这个使用函数 sinc 的程序，可用来计算上述近似滤波器的冲激响应。 clf; fc = 0.25; n = -6.5:1:6.5; y = 2*fc*sinc(2*fc*pi*n);k = n+6.5; 第 19 页 共 41 页 stem(k,y);title(N = 13);axis(0 13 -0.2 0.6); xlabel(Time index n);ylabel(Amplitude);grid; 实验内容 42 运行例 41 的程序计算并画出近

36、似理想低通滤波器的冲激响应。 在例 41 中的低通有限冲激响应滤波器的长度是多少？该程序中哪个语句 确定滤波器的长度？哪个参数控制截止频率？ 修改例 41 中的程序，计算长度为 20，截止角频率为 0.45 的有限冲激响应 低通滤波器的冲激响应。 （hn不变） 3. 2例 43 ：线性相位有限冲激响应滤波器时不变系统的零点的位置研究。 下列程序先产生四类滤波器类线性时不变系统（有限长线性相位滤波器）的冲 激响应，然后产生零极点图，最终显示零点的位置。 clf; b = 1 -8.5 30.5 -63; num1 = b 81 fliplr(b); num2 = b 81 81 fliplr(b

37、); num3 = b 0 -fliplr(b); num4 = b 81 -81 -fliplr(b); n1 = 0:length(num1)-1; n2 = 0:length(num2)-1; subplot(2,2,1); stem(n1,num1); xlabel(Time index n);ylabel(Amplitude); grid; title(Type 1 FIR Filter); subplot(2,2,2); stem(n2,num2); xlabel(Time index n);ylabel(Amplitude); grid; title(Type 2 FIR Fil

38、ter); subplot(2,2,3); stem(n1,num3); xlabel(Time index n);ylabel(Amplitude); grid; title(Type 3 FIR Filter); subplot(2,2,4); stem(n2,num4); xlabel(Time index n);ylabel(Amplitude); grid; title(Type 4 FIR Filter); pause subplot(2,2,1); zplane(num1,1); title(Type 1 FIR Filter); subplot(2,2,2); zplane(n

39、um2,1); title(Type 2 FIR Filter); subplot(2,2,3); zplane(num3,1); title(Type 3 FIR Filter); subplot(2,2,4); zplane(num4,1); title(Type 4 FIR Filter); disp(Zeros of Type 1 FIR Filter are); disp(roots(num1); disp(Zeros of Type 2 FIR Filter are); disp(roots(num2); disp(Zeros of Type 3 FIR Filter are);

40、第 20 页 共 41 页 disp(roots(num3); disp(Zeros of Type 4 FIR Filter are); disp(roots(num4); 4实验问题回答 （1）运行例 4-1 的程序，画出给定系统频率响应的幅度和相位响应。 （2）运行例 4-2 的程序，画出这个滤波器的幅度和相位响应，并验证它是一 个低通滤波器。 （3）运行例 43 的程序，生成每一类线性相位有限冲激响应滤波器的冲激响 应。并回答每一个有限冲激响应滤波器的长度是多少？ （4）在例 43 中验证冲激响应序列的对称性，每一个滤波器的零点位置及它 们的线性相位特性。并回答这些滤波器的群延迟是多少

41、？ 四实验报告要求 1按照实验内容要求完成相关实验程序，并得出相关的实验结果（包括图 形结果） 。 2回答实验中提出的问题。 3总结本次实验结果，按照实验报告格式要求，书写实验报告。 五实验设备 PC 机，MATLAB 软件 第 21 页 共 41 页 实验五：实验五： IIR 数字滤波器的设计数字滤波器的设计 一实验目的 1加深对 IIR 滤波器的理解 2学会用直接法和间接法设计满足响应要求的 IIR 数字滤波器 二实验相关知识准备 1模拟滤波器的设计函数 为了利用模拟滤波器设计数字滤波器，先必须设计对应的模拟滤波器，常用的 模拟滤波器有：bessel 滤波器，butterworth 滤波器

42、，chebyshevI 型滤波器， chebyshevII 型滤波器及椭圆函数滤波器。 1 模拟滤波器的设计函数 第 22 页 共 41 页 设计 bessel 模拟低通滤波器 z,p,k=besselap(n) %设计 bessel 模拟低通滤波器 2） 设计 butterworth 模拟低通滤波器 z,p,k=buttap(n) %设计 butterworth 模拟低通滤波器 3） 设计 chebyshevI 型模拟低通滤波器 z,p,k=cheb1ap(n,Rp) %设计 chebyshevI 型模拟低通滤波器 %Rp：通带内的波纹系数，单位分贝 4） 设计 chebyshevII 型模

43、拟低通滤波器 z,p,k=cheb2ap(n,Rs) %设计 chebyshevII 型模拟低通滤波器 %Rs：阻带内的波纹系数低于通带 Rs 分贝 5）设计椭圆模拟滤波器 z,p,k=ellipap(n,Rp.Rs) %设计椭圆模拟滤波器 2滤波器阶数的选择 下列函数除了能选择模拟滤波器的阶数外，同时也能选择数字滤波器的阶数。 选择 butterworth 滤波器阶数 数字域：n,Wn=buttord(Wp,Ws,Rp,Rs) 模拟域：n,Wn=buttord(Wp,Ws,Rp,Rs,s) 2）选择 chebyshevI 型滤波器阶数 数字域：n,Wn=cheb1ord(Wp,Ws,Rp,R

44、s) 模拟域：n,Wn=cheb1ord(Wp,Ws,Rp,Rs,s) 3）选择 chebyshevII 型滤波器阶数 数字域：n,Wn=cheb2ord(Wp,Ws,Rp,Rs) 模拟域：n,Wn=cheb2ord(Wp,Ws,Rp,Rs,s) 选择椭圆滤波器阶数 数字域：n,Wn=ellipord(Wp,Ws,Rp,Rs) 模拟域：n,Wn=ellipord(Wp,Ws,Rp,Rs,s) 注意： n：返回符合要求性能指标的数字滤波器或模拟滤波器的最小阶数 Wn：滤波器的截至频率（即 3db 频率） Wp：通带的截至频率，Ws：阻带的截至频率，单位 rad/s。且均为规一化频 第 23 页

45、共 41 页 率，即。1 对应 弧度。1)(0WsWp 频率规一化：信号处理工具箱中使用的频率为莱奎斯特频率，根据香农定理， 它为采样频率的一半，在滤波器设计中的截止频率均使用莱奎斯特频率进行规一化。 规一化频率转换为角频率，则将规一化频率乘以 。如果将规一化频率转换为 Hz，则将规一化频率乘以采样频率的一半。 3 模拟频率变换 1）低通到低通模拟滤波器变换 bt,at=lp2lp(b,a,Wo) lp2lp 函数将截至频率为 1rad/s（规一化截至频率）的模拟低通滤波器原型变换 成截至频率为 Wo 的低通滤波器。 2）低通到带通模拟滤波器变换 bt,at=lp2bp(b,a,Wo,Bw)

46、lp2bp 函数可将截止频率为 1rad/s 的模拟低通滤波器原型变换成具有指定带宽 Bw 和中心频率 Wo 的带通滤波器。其中 Wo=sqrt(W1*W2)，Bw=W2-W1 ，W2 高端截止频率，W1 低端截止频率。 3）低通到高通模拟滤波器变换 bt,at=lp2hp(b,a,Wo) 4）低通到带阻模拟滤波器变换 bt,at=lp2bs(b,a,Wo,Bw) 5）besself 模拟滤波器设计 b,a=besself(n,Wn) b,a=besself(n,Wn,ftype) z,p,k=besself(n,Wn) z,p,k=besself(n,Wn,ftype) 2由模拟滤波器变换成

47、等效的数字滤波器 方法一：双线性变换法实现模拟滤波器到数字滤波器的变换 双线性变换将 S 域映射成 Z 域，从而将模拟滤波器变换成等效的数字滤波器 zd,pd,kd=bilinear(z,p,k,Fs)为零极点增益表示的 bilinear 函数，其中 z,p,k 为 S 域传 递函数的零点、极点、和增益，Fs 为取样频率，zd,pd,kd 为双线性变换后 Z 域传递 函数的零点、极点和增益。 第 24 页 共 41 页 numd,dend=bilinear(num,den,Fs,Fp)为传输函数表示的 bilinear 函数 方法二：冲激响应不变法实现模拟到数字滤波器的变换 bz,az=imp

48、invar(b,a,Fs) bz,az=impinvar(b,a)采用默认值为 1hz 的 Fs。 小结：小结： 设计 IIR 数字滤波器的一般步骤： 把给出的数字滤波器的性能指标转换为模拟滤波器的性能指标 根据转换后的性能指标，通过滤波器阶数选择函数，来确定滤波器的最小阶 数 N 和固有频率 Wn 由最小阶数 N 得到低通滤波器原型 由固有频率 Wn 把模拟低通滤波器原型转换为低通、高通、带通、带阻滤波 器 运用脉冲响应不变法或双线性变换法把模拟滤波器转换成数字滤波器 MATLAB 信号处理工具箱提供了几个用于直接设计 IIR 数字滤波器的函数， 这些函数把上面这些复杂的步骤融合为一个整体，

49、为设计滤波器带来了极大的方便。 3直接设计 IIR 数字滤波器 1Butterworth 模拟和数字滤波器设计 数字域：b,a=butter(n,Wn)可设计出截止频率为 Wn 的 n 阶 butterworth 滤波器 b,a=butter(n,Wn,ftype)当 ftype=high 时，可设计出截止频率为 Wn 的高 通滤波器；当 ftypestop 时，可设计出带阻滤波器 z,p,k=butter(n,Wn) zp,k=buter(n,Wn,ftype) A,B,C,D=butter(n,Wn) A,B,C,D=butter(n,Wn,ftype) 模拟域：b,a=butter(n,

50、Wn,s)可设计出截止频率为 Wn 的 n 阶模拟 butterworth 滤 波器， 其余形式类似于数字域的。 2chebyshevI 型滤波器（通带等波纹）设计 数字域：b,a=cheby1(n,Rp,Wn)可设计出 n 阶 chebyshevI 滤波器,其截止频率由 Wn 第 25 页 共 41 页 确定，通带内的波纹由 Rp 确定 b,a=cheby1(n,Rp,Wn,ftype)当 ftype=high 时，可设计出截止频率为 Wn 的高通滤波器；当 ftypestop 时，可设计出带阻滤波器 z,p,k=cheby1(n,Rp,Wn) zp,k= cheby1 (n,Rp,Wn,f

51、type) A,B,C,D= cheby1 (n,Rp,Wn) A,B,C,D= cheby1 (n,Rp,Wn,ftype) 模拟域：b,a= cheby1 (n,Rp,Wn,s)可设计出截止频率为 Wn 的 n 阶 chebyshevI 型 模拟滤波器，其余形式类似于数字域的。 3chebyshevII 型滤波器（阻带等波纹）设计 数字域：b,a=cheby2(n,Rs,Wn)可设计出 n 阶 chebyshevI 滤波器,其截止频率由 Wn 确定，阻带内的波纹由 Rs 确定 b,a=cheby2(n,Rs,Wn,ftype)当 ftype=high 时，可设计出截止频率为 Wn 的高通滤

52、波器；当 ftypestop 时，可设计出带阻滤波器 z,p,k=cheby2(n,Rs,Wn) zp,k= cheby2(n,Rs,Wn,ftype) A,B,C,D= cheby2 (n,Rs,Wn) A,B,C,D= cheby2(n,Rs,Wn,ftype) 模拟域：b,a= cheby2(n,Rs,Wn,s)可设计出截止频率为 Wn 的 n 阶 chebyshevII 型模拟滤波器，其余形式类似于数字域的。 4椭圆滤波器设计 数字域：b,a=ellip(n,Rp,Rs,Wn)可设计出 n 阶椭圆滤波器,其截止频率由 Wn 确定， 通带内的波纹由 Rp 确定,阻带内的波纹由 Rs 确定

53、。 b,a=ellip(n,Rp,Rs,Wn,ftype)当 ftype=high 时，可设计出截止频率为 Wn 的高通滤波器；当 ftypestop 时，可设计出带阻滤波器 z,p,k=ellip(n,Rp,Rs,Wn) zp,k=ellip(n,Rp,Rs,Wn,ftype) A,B,C,D=ellip(n,Rp,Rs,Wn) A,B,C,D=ellip(n,Rp,Rs,Wn,ftype) 模拟域：b,a=ellip(n,Rp,Rs,Wn,s)可设计出截止频率为 Wn 的 n 阶椭圆型模拟 滤波器，其余形式类似于数字域的。 第 26 页 共 41 页 54基于幅度特性设计 IIR 数字滤波

54、器 b,a=yulewalk(n,f,m)可得到以 b,a 表示的 n 阶 IIR 滤波器，其幅频特性逼近由频率 矢量 f 和幅度矢量 m 所给定的特性。矢量 f 和 n 的长度必须相同。 65基于频率特性设计 IIR 数字滤波器 b,a=invfreqz(h,w,nb,na)可得到以 b,a 表示的 IIR 滤波器，其幅频特性逼近由频率矢 量 w 和复频率响应矢量 h 所给定的特性。矢量 w 和 h 的长度必须相同，nb,na 分别 为指定的分子和分母多项式的阶数。 三实验内容 1 例 51 1设计一 butterworth 低通滤波器，通带范围 0100hz，通带波纹小于 3db，阻带 为

55、-30db，并利用最小的阶数来实现。采样频率为 1000hz。 wp=100/500;ws=200/500; %频率归一化,ws? n,wn=buttord(wp,ws,3,30); b,a=butter(n,wn); freqz(b,a,512,1000) 例 5-2 用双线性变换法设计一个切比雪夫 I 型 IIR 数字高通滤波器 高通滤波器的指标为：通带截止频率 Fp=700HZ, 阻带截止频率 Fs=500HZ, 通带波纹 p=1dB,最小阻带衰减 s=32Db,抽样频率为 Fr=2k HZ. 解：(1)归一化截止频率为： p=0.7 s=0.5 (2)双线性变换的频率预畸变(模拟高通滤

56、波器) p=1.9626105 s=1.0 (3) 归一化频率变换（模拟高通-模拟低通）要求与 1相同，但采用 椭圆函数滤波器实现。 wp=100/500;ws=200/500;%频率归一化 n,wn=ellipord(wp,ws,3,30); b,a=ellip(n,3,30,wn); 第 27 页 共 41 页 freqz(b,a,512,1000) 应用最小二乘拟合的方法设计一个 8 阶的低通滤波器，指定的频率点与 对应的幅度分别由 f 和 m 给出，并画出它希望的频率响应和实际的频率响应。 f=0 0.6 0.6 1; m=1 1 0 0; b,a=yulewalk(8,f,m); h

57、,w=freqz(b,a,128); plot(f,m,w/pi,abs(h),-); legend(ideal,yulewalk designed); title(comparison of frequency response magnitudes) 选择 p=1 则 s=1.9626105 (4) 设计相应的模拟低通滤波器（以下步骤参看下面程序） (5) 模拟低通滤波器到模拟高通滤波器转换 (6) 双线性变换完成模拟高通滤波器转换成数字高通滤波器. 程序 5-2： clf; N,Wn=cheb1ord(1,1.9626105,1,32,s); B,A = cheby1(N,1,Wn,s)

58、; BT,AT =lp2hp(B,A, 1.9626105); num,den =bilinear(BT,AT,0.5); 2实验内容 53 确定通带截止频率为 500HZ，阻带截止频率为 600HZ 以及通带的最大衰减量 为 1db，阻带的最小衰减量为 50db，采样频率为 2000HZ 的数字低通巴特沃斯滤波 器的阶数，并画出其幅频特性图。 实验内容 54 采用直接设计法设计一个通带截止频率为 500HZ800HZ、阻带截止频率为 600HZ700HZ、通带内波纹的最大衰减为 0.1db、阻带内波纹的最小衰减为 50db、采样频率为 2000HZ 的带阻椭圆数字滤波器，并画出其幅频响应图。

59、 第 28 页 共 41 页 实验内容 55 模拟信号 xa（t）=sin(500*pi*t)+cos(50*pi*t) 有系统完成xa（t） A/D H(z) D/A ya(t)的处理， 采样频率为 1000HZ。 采用双线性变换法设计一个高通切比雪夫 2 型数字滤波器，要求所设计的滤波器阶 数最小，以小于 1db 的衰减通过 150hz 的分量，以至少 40db 抑制 100hz 分量。滤 波器应有单调的通带和等波纹的阻带。求此滤波器并画出幅频响应图（单位 db） 产生上述信号xa（t）的 150 个样本点，用图形文件在时域和频域上分别表示该信 号。 然后通过上述所设计的滤波器对（2）所产

60、生的信号滤波，并分别在时域和频域上 结合滤波后的输出信号进行结果分析。 3实验问题回答 （1）完成例 51、52 的程序并观察 IIR 滤波器的特点。 （2）在实验内容 5-3、5-4 和 55 中按要求完成程序，并回答提出的问题。 四实验报告要求 1按照实验内容要求完成相关实验程序，并得出相关的实验结果（包括图 形结果） 。 2回答实验中提出的问题。 3总结本次实验结果，按照实验报告格式要求，书写实验报告。 五实验设备 PC 机，MATLAB 软件 第 29 页 共 41 页 实验六：实验六： FIR 数字滤波器的设计数字滤波器的设计 一实验目的 1加深对 FIR 滤波器的理解 2学会用窗函