数字信号处理实验指导书.doc

数字信号处理实验指导书

———————————————————————————————— 作者：
———————————————————————————————— 日期：

2

个人收集整理 勿做商业用途

数字信号处理
实 验 指 导 书

武汉理工大学教材中心
2012年7月

45

实验一 时域离散信号的产生

一、实验目的
1、了解常用时域离散信号及其特点;
2、掌握MATLAB程序的编程方法；
3、熟悉MATLAB函数的调用方法。
二、实验原理
在时间轴上的离散点取值的信号，称为离散时间信号.离散时间信号只在某些离散的瞬时给出函数值，而在其他时刻无定义。它是时间上不连续按一定先后次序排列的一组数的集合，称为时间序列，用x（n)表示，n取整数代表时间的离散时刻。
在MATLAB中用向量来表示一个有限长度的序列。
常用离散信号：
1、单位抽样序列

2、单位阶跃序列

3、实指数序列

4、复指数序列

5、正(余）弦序列

6、随机序列
在利用计算机进行系统的研究时，经常需要产生随机信号,MATLAB提供一个工具函数rand来产生随机信号。
7、周期序列

三、实验用函数
1、stem
功能：绘制二维图形.
调用格式：
stem（n,x）；n为横轴，x为纵轴的线性图形。
2、length
功能：计算某一变量的长度或采样点数。
调用格式：
N=length（t);计算时间向量t的个数并赋给变量N。
3、axis
功能：限定图形坐标的范围。
调用格式：
axis(［x1,x2，y1,y2］);横坐标从x1-x2，纵坐标从y1—y2。
4、zeros
功能：产生一个全0序列。
调用格式：
x=zeros（1,n)；产生n个0的序列.
5、ones
功能：产生一个全1序列。
调用格式：
y=ones（1，n）；产生n个1的序列.
四、参考实例
例1.1 用Matlab产生单位抽样序列。
%先建立函数impseq（n1，n2，n0）
function [x，n]=impseq（n1,n2,n0)
n=[n1:n2];
x=［(n—n0)==0];
%编写主程序调用该函数
[x，n］=impseq（—2，8，2);
stem(n,x）
程序运行结果如图1-1所示:

图1—1 单位抽样序列
例1.2实数指数序列（运算符“.^”）
Matlab程序如下:
n=[0:10］;
x=0。9.^n；
stem（n，x)
程序运行结果如图1-2所示

图1—2 实数指数序列
例1。3复数指数序列（）
Matlab程序如下：
n=[-10:10]; alpha=-0。1+0.3＊j； x=exp（alpha*n)；
real_x=real（x)； image_x=imag（x);
mag_x=abs（x)； phase_x=angle(x）；
subplot（2，2,1)； stem（n，real_x)
subplot（2，2，2）； stem(n,image_x)
subplot（2，2，3); stem（n,mag_x）
subplot（2,2，4）; stem（n，phase_x）
程序运行结果如图1—3所示

图1—3 复数指数序列
例1。4正、余弦序列（）
Matlab程序如下:
n=［0：10]；
x=3＊cos（0.1＊pi＊n+pi/3);
stem(n，x)
程序运行结果如图1—4所示

图1—4 正、余弦序列
例1。5随机序列
rand（1,N)产生其元素在[0，1］之间均匀分布长度为N的随机序列
randn(1，N)产生均值为0，方差为1，长度为N的高斯随机序列
例1.6周期序列
如何生成周期序列
1、 将一个周期复制p次;
2、借助矩阵运算、matlab下标能力。先生成一个包含p列x（n）值的矩阵，然后用结构（:）来把p列串接成一个长周期序列。因为这个结构只能用于列向,最后还需要做矩阵转置获得所需序列。
Matlab程序如下：
x=[1,2,3]； ％一个x（n)
xn=x’＊ones(1,3) %生成p列x(n）
xn=xn（：)' %将p列串接成长列序列并转置
stem(xn)
程序运行的结果如图1-5所示

图1—5 周期序列
五、实验任务
1、调试部分例题程序，掌握Matlab基本操作方法。
2、编写程序，完成下列函数波形：
1）利用zeros函数生成单位抽样序列;
2）利用zeros函数和ones函数生成单位阶跃序列；

六、实验报告
1、简述实验目的、原理。
2、写出上机调试通过的实验任务的程序并描述其图形曲线.

实验二 离散序列的基本运算
一、实验目的
1、加强MATLAB运用。
2、了解离散时间序列在时域中的基本运算。
3、熟悉相关函数的使用方法，掌握离散序列运算程序的编写方法。
二、实验原理
离散序列的时域运算包括信号的相加、相乘,信号的时域变换包括信号的移位、反折、倒相及尺度变换等.
在MATLAB中，序列的相加和相乘运算是两个向量之间的运算，因此参加运算的两个序列必须具有相同的长度，否则不能直接进行运算，需要进行相应的处理后再进行运算。三、实验用函数
1、find
功能：寻找非零元素的索引号。
调用格式：
find(（n〉=min(n1））＆（n〈=max（n1)))：在符号关系运算条件的范围内寻找非零元素的索引号。
2、fliplr
功能：对矩阵行元素进行左右翻转。
调用格式:
x1=fliplr(x)：将x的行元素进行左右翻转,赋给变量x1。
四、实例
1、信号的时域变换
1）序列移位
将一个离散序列进行移位,形成新的序列：x1(n）=x（n-m)。当m〉0时，原序列向右移m位，当m〈0时，原序列向左移.
%建立移位函数(sigshift（x，m,n0)）
function [y，n]=sigshift（x,m，n0)
n=m+n0；
y=x；
2）序列反折
在这个运算中，x（n)以n=0为基准点，以纵轴为对称轴反折得到一个新的序列。
y(n)=｜x(-n）|
在MATLAB中提供了fliplr函数实现序列反折。
%建立反折函数(sigfold(x，n)）
function ［y，n]=sigfold(x,n）
y=fliplr（x);
n=-fliplr（n)；
3）序列倒相
是求一个与原序列的向量值相反,对应的时间向量不变的新序列。
4)序列的尺度变换
通过对时间轴的放大或压缩形成新的序列.
2、序列的算术运算
1）序列相加
序列相加是指两个序列中相同序号的序列值逐项对应相加，形成新的序列。
参加运算的两个序列的维数不同时
％建立通用函数
function ［y,n]=sigadd(x1,n1，x2，n2)
n=min（min（n1),min(n2)) : max（max(n1），max(n2)）；
y1=zeros（1，length（n))；
y2=y1;
y1（find((n>=min（n1）)＆（n〈=max(n1））==1))=x1;
y2（find（（n>=min（n2))&(n<=max(n2）)==1））=x2；
y=y1+y2；
1）序列相乘
序列相加是指两个序列中相同序号的序列值逐项对应相乘,形成新的序列。
参加运算的两个序列的维数不同时处理方法与序列相加相同。

五、实验任务
1、理解序列运算的性质，了解函数语句的意义。
2、利用例题函数完成下列序列运算
1）已知x1(n）=u（n+1） （-3<n〈5);
x2（n)=u(n—3） (-4〈n<7）
求:x（n）=x1(n)+x2(n)
2）已知x1(n)=3e-0。25n (—2〈n〈8)
x2（n)=u(n+1) （-3〈n<6)
求：x(n）=x1(n)＊x2（n）
六、实验报告
1、简述实验目的和原理。
2、列写上机调试通过的程序，并描绘其波形曲线。

实验三 离散卷积的原理及应用
一、实验目的
1、通过实验进一步理解卷积定理,了解卷积过程；
2、掌握应用线性卷积求解离散时间系统响应的基本方法。
二、实验原理
对于线性移不变离散系统，任意的输入信号x（n)可以用及其位移的线性组合来表示，即

当输入为时，系统的输出y(n)=h(n),由系统的线性移不变性质可以得到系统对x（n）的响应y(n）为

称为离散系统的线性卷积，简写为
y（n)=x(n)＊h（n）
也就是说，如果已知系统的冲激响应,将输入信号与系统的冲激响应进行卷积运算，即可求得系统的响应。
三、实验用函数
1、卷积函数conv
功能：进行两个序列的卷积运算.
调用格式：
y=conv（x，h)；用于求解两有限长序列的卷积。
2、sum
功能：求各元素之和。
调用格式：
y=sum（x);求序列x中各元素之和。
3、hold
功能：控制当前图形窗口是否刷新的双向切换开关。
调用格式:
hold on:使当前图形窗口中的图形保持且不被刷新，准备接受绘制新的图形。
hold off:使当前图形窗口中的图形不具备不被刷新的性质。
4、impz
功能：求解数字系统的冲激响应。
调用格式:
［h，t]=impz(b，a）;求解数字系统的冲激响应h，取样点数为缺省值.
[h,t]=impz（b,a,n）；求解数字系统的冲激响应h，取样点数由n确定.
impz(b,a)；在当前窗口用stem（t，h）函数绘制图形。
5、dstep
功能：求解数字系统的阶跃响应。
调用格式:
[h,t］=dstep（b，a）；求解数字系统的冲激响应h，取样点数为缺省值。
［h，t］=dstep（b，a,n)；求解数字系统的冲激响应h,取样点数由n确定。
dstep（b，a）；在当前窗口用stairs（t，h)函数绘制图形。
四、参考实例
在利用Matlab提供的卷积函数进行卷积运算时，主要是确定卷积结果的时间区间。conv函数默认两信号的时间序列从n=0开始，卷积结果对应的时间序列也从n=0开始。如果信号不是从0开始，则编程时必须用两个数组确定一个信号，其中，一个数组是信号波形的幅度样值，另一个数组是其对应的时间向量。
％建立一个适用于信号从任意时间开始的通用函数
function [y,ny]=sconv（x，h,nx,nh,p)
y=conv(x，h）；
n1=nx（1)+nh(1); %计算y的非零样值的起点位置
n2=nx(length（x）)+nh(length(h））； %计算y的非零样值的宽度
ny=n1:p:n2； ％确定y的非零样值的时间向量
五、实验任务
已知一个IIR数字低通滤波器的系统函数公式为

输入一个矩形信号序列 x=square(n/5) (-2〈n〈10）,求该系统的响应。
六、实验报告
1、简述实验目的、原理。
2、写出上机调试通过的实验任务的程序并描述其图形曲线。

实验四 离散傅立叶级数
一、实验目的
1、加深对离散周期序列傅里叶级数基本概念的理解。
2、掌握MATLAB求解周期序列傅里叶级数变换和逆变换的方法。
3、观察离散周期序列的重复周期数对频谱特性的影响。
二、实验原理
离散时间序列x(n)满足x(n)=x(n+rN），称为离散周期序列,其中N为周期，x(n）为主值序列.
周期序列可用离散傅里叶级数表示成
n=0，1,…,N—1
其中，是周期序列离散傅里叶级数第K次谐波分量的系数,也称为周期序列的频谱，可表示为
k=0,1，…，N-1
上面两式是周期序列的一对傅里叶级数变换对.
令,以上两式可简写为：

三、实验用函数
1、mod
功能:模除求余。
调用格式:
mod（x,m）：x整除m取正余数。
2、floor
功能：向—舍入为整数。
调用格式：
floor(x）：将x向-舍入为整数。
四、实例
1、周期序列的傅里叶变换和逆变换
依据变换公式编写通用函数
1）离散傅里叶级数正变换通用函数
function xk=dfs（xn，N)
n=[0：1:N—1］； ％n的行向量
k=n； %k的行向量
WN=exp(-j*2*pi/N)； %WN因子
nk=n’*k; %产生一个含nk值的N乘N维矩阵
WNnk=WN.^nk； %DFS矩阵
xk=xn* WNnk； %DFS系数行向量
2）离散傅里叶级数逆变换通用函数
function xn=idfs(xk，N)
n=［0：1：N—1］; %n的行向量
k=n; %k的行向量
WN=exp(—j＊2＊pi/N)； %WN因子
nk=n'*k； ％产生一个含nk值的N乘N维矩阵
WNnk=WN.^（—nk); ％DFS矩阵
xn=xk* WNnk/N； ％DFS系数行向量
例:已知一个周期性矩形序列的脉冲宽度占整个周期的1/4，一个周期的采样点为16点。用傅里叶级数求信号的幅度和相位频谱;求傅里叶级数逆变换的图形，与原信号图形进行比较。
MATLAB程序
N=16;
xn=［ones（1,N/4）,zeros（1，3＊N/4)];
n=0:N-1;
xk=dfs（xn，N）;
xn1=idfs（xk，N)；
subplot(2，2，1）; stem(n,xn）; title(’x（n)');
subplot(2,2，2); stem（n,abs（xn1）); title('idfs(|X（k)|)’）;
subplot(2，2,3)； stem(n,abs（xk)）; title（'|X(k）｜’）;
subplot(2，2，4); stem(n，angle（xk)); title（'arg｜X（k）|’);
程序运行结果如图4-1

图4—1
2、周期重复次数对序列频谱的影响
理论上讲，周期序列不满足绝对可积条件，要对周期序列进行分析，可以先取K个周期进行处理,然后让K无限增大,研究其极限情况。这样可以观察信号序列由非周期到周期变换时,频谱由连续谱逐渐向离散谱过渡的过程。
例：已知一个矩形序列的脉冲宽度占整个周期的1/2，一个周期的采样点数为10，用傅立叶级数变换求信号的重复周期数分别为1、4、7、10时的幅度频谱.
MATLAB程序：
xn=［ones(1，5),zeros（1,5）]；
Nx=length（xn）；
Nw=1000;dw=2＊pi/Nw;
k=floor(（-Nw/2+0.5)：（Nw/2+0.5)）;
for r=0：3；
K=3＊r+1;
nx=0：(K*Nx—1)；
x=xn（mod(nx，Nx）+1）；
Xk=x*（exp（—j*dw*nx'＊k)）/K;
subplot(4,2，2*r+1）; stem(nx,x）
axis（［0,K＊Nx-1,0,1.1］）； ylabel（'x（n)'）;
subplot(4,2,2＊r+2)； plot（k＊dw,abs（Xk)）
axis([—4,4,0,1.1*max（abs(Xk))]); ylabel（'X(k）'）;
end
程序运行结果如图4—2

图4-2
从上图可以看出，信号序列的周期数越多，则频谱越是向几个频点集中，当信号周期数趋于无穷大时 ，频谱转化为离散谱。
五、实验任务
1、输入并运行例题程序，熟悉基本指令的使用。
2、已知一个信号序列的主值为x(n)=［0,1，2，3,2,1，0］,显示两个周期的信号序列波形.
求解：用傅里叶级数求信号的幅度和相位频谱;求傅里叶级数逆变换的图形，与原信号图形进行比较。
六、实验报告
1、简述实验目的、原理.
2、写出上机调试通过的实验任务的程序并描述其图形曲线。

实验五 离散傅里叶变换
一、实验目的
1、加深对离散傅里叶变换基本概念的理解。
2、了解有限长序列傅里叶变换与周期序列傅里叶级数的联系。
3、熟悉相关函数的使用方法。
二、实验原理
有限长序列的傅里叶变换和逆变换
对于非周期序列,在实际中常常使用有限长序列。有限长序列x（n）表示为

x（n)是非周期序列，但可以理解为某一周期序列的主值序列。由离散傅立叶级数DFS和IDFS引出有限长序列的离散傅立叶正、逆变换关系式.

DFT与DFS的关系
比较两者的变换对，可以看出两者的区别仅仅是将周期序列换成了有限长序列.
有限长序列x(n)可以看作是周期序列的一个周期；反之周期序列可以看作是有限长序列x(n)以N为周期的周期延拓.
由于公式非常相似,在程序编写上也基本一致.
三、实例
1、已知有限长序列x（n)为:
x(n）=［0,1,2,3,4，5,6，7，8,9]，求x(n）的DFT和IDFT。要求
1）画出序列傅里叶变换对应的｜X(k)｜和arg［X(k）]图形.
2)画出原信号与傅里叶逆变换IDFT[X（k）］图形进行比较。

MATLAB程序：
xn=［0，1，2，3，4,5，6,7，8,9］； N=length(xn)；
n=0:N-1; k=0：N—1;
xk=xn＊exp（—j＊2*pi/N）.^（n’＊k）;
x=(xk＊exp（j*2＊pi/N）.^（n'*k)）/N；
subplot（2，2，1）; stem（n，xn); title(‘x(n）’)；
subplot(2,2，2）； stem(n，abs（x））; title（‘IDFT|X(k）|')；
subplot(2，2，3); stem（k,abs（xk））; title（‘｜X（k)｜’）；
subplot(2，2，4）; stem（k，angle（Xk)）; title(‘arg|X(k)｜')；
程序运行结果如图5-1：

图5-1
由上图可看出，与周期序列不同的是，有限长序列本身是仅有N点的离散序列,相当于周期序列的主值部分。因此,其频谱也对应序列的主值部分，是含N点的离散序列。
2、有限长序列DFT与周期序列DFS的联系
已知周期序列的主值x（n)=［0,1，2,3,4，5］，求x（n)周期重复次数为4次时的DFS.要求
1)画出原主值序列和信号周期序列；
2）画出序列傅里叶变换对的图形.
MATLAB程序：
xn=［0,1，2，3,4,5］; N=length（xn）；
n=0：4*N-1； k=0:4＊N-1;
xn1=xn(mod（n,N)+1）;
xk=xn1*exp（—j*2*pi/N)。^（n’＊k）；
subplot(2,2,1)； stem(xn)； title（’原主值信号x（n）')；
subplot（2，2，2）; stem(n，xn1）； title（'周期序列信号'）;
subplot(2,2，3）; stem(k，abs（xk））; title('｜X(k）｜’);
subplot（2，2,4）； stem（k，angle（xk)）； title（'arg｜X（k）｜')；
程序运行结果如图5-2：

图5—2
与上一个例题比较,有限长序列x（n）可以看成是周期序列的一个周期,反之，周期序列可以看成是有限长序列以N为周期的周期延拓。频域上的情况也是相同的。从这个意义上说,周期序列只是有限个序列值有意义.
四、实验任务
1、输入并运行例题程序,熟悉基本指令的使用.
2、验证离散傅里叶变换的线性性质。
有两个有限长序列分别为x1（n）和x2(n)，长度分别为N1和N2，且y（n)=ax1(n）+bx2(n)，（a,b均为常数）,则该y（n）的N点DFT为
Y（k)=DFT［y(n)］=aX1（k）+bX2（k） （0〈=k<=N—1）
其中：N=max(N1,N2)，X1（k）和X2(k)分别为x1(n）和x2（n）的N点DFT。
已知序列: x1（n）=［0,1，2，4]
x2(n）=[1，0，1,0,1]
五、实验报告
1、简述实验目的、原理。
2、写出上机调试通过的实验任务的程序并描述其图形曲线.

实验六 快速傅里叶变换
一、实验目的
1、加深对快速傅里叶变换基本理论的理解。
2、了解用MATLAB语言进行快速傅里叶变换的方法。
3、掌握常用函数的调用方法.
二、实验原理
DFT是在时域和频域均为离散序列的变换方法，它适用于有限长序列。但如果按照变换公式进行运算的话，当序列长度很大时，将占用很大的内存空间，运算时间也会很长,无法实时处理问题。
快速傅里叶变换是用于提高DFT运算的高速运算方法的统称,FFT是其中的一种，FFT不是一种新的变换形式，它仅仅只是一种快速算法。FFT主要有时域抽取算法和频域抽取算法,基本思想是将一个长度为N的序列分解成多个短序列再进行运算，如基2算法、基4算法等等，从而可以大大缩短运算时间。
三、实验用函数
1、fft
功能：一维基2快速傅里叶变换.
调用格式:
y=fft(x)：利用FFT算法计算矢量x的离散傅里叶变换。当x为2的幂次方时，采用高速基2FFT算法，否则为稍慢的混合算法。
y=fft（x,N)：采用n点FFT。当x的长度小于N时，FFT函数在x的尾部补零，以构成N点数据,当x的长度大于N时，FFT函数在x的尾部截断，以构成N点数据。
2、ifft
功能：一维基2快速傅里叶逆变换。
调用格式:
与FFT调用方法相同,只需改换函数名.
3、fftshift
功能：对FFT的输出进行重新排列，将零频分量移到频谱的中心.
调用格式：
y=fftshift（x）:对FFT的输出进行重新排列，将零频分量移到频谱的中心。
四、实例
1、用MATLAB工具箱函数FFT进行频谱分析时需要注意：
1）函数fft的返回值y的数据结构的对称性
若已知序列x=［4,3，2，6,7，8，9,0]，求X（k）=DFT[x（n）］
利用函数fft计算,其MATLAB程序如下:
N=8;
n=0：N-1；
xn=[4,3，2，6，7,8，9,0］;
xk=fft(xn)'
程序运行结果如下:
xk =
39.0000
—10。7782 - 6。2929i
0 + 5.0000i
4。7782 + 7。7071i
5.0000
4。7782 — 7.7071i
0 — 5。0000i
—10。7782 + 6。2929i
由程序运行结果可见,xk的第一行元素对应频率值为0,第五行元素对应频率为莱奎斯特（Nyquist）频率，即标准频率值为1。因此第一行至第五行对应的标准频率为0~1。而第五行至第八行对应的是负频率,其K（x）值是以Nyquist频率为轴对称.
一般情况，对于N点的x(n)序列的FFT是N点的复数序列，其点n=N/2+1对应Nyquist频率,作谱分析时仅取序列X(k）的前一半即可，其后一半序列和前一半是对称的。
2）频率计算
若N点序列x（n)（n=0,1，…，N-1)是在采样频率fs(Hz）下获得的。它的FFT也是N点序列，即X（k）（k=0，1，…，N-1),则第K点对应实际频率值为：
f=k＊fs/N
3）作FFT分析时，幅值大小与FFT选择的点数有关，但不影响分析结果。
2、已知信号由15Hz幅值0.5的正弦信号和40Hz幅值2的正弦信号组成，数据采样频率为100Hz，试绘制N=128点DFT的幅频图。
MATLAB程序如下：
fs=100;
N=128;
n=0：N—1;
t=n/fs；
x=0。5＊sin（2＊pi*15*t）+2＊sin(2*pi*40*t)；
y=fft（x,N）;
f=(0:length（y）—1）’*fs/length（y）；
mag=abs(y）；
stem(f,mag);
title（'N=128点')
程序运行的结果如图6-1

图6-1
如图所示，由于信号采样频率为100Hz，故其莱奎斯特频率为50Hz，图中整个频谱图是以莱奎斯特频率为轴对称的.因此利用FFT对信号作谱分析时，只要考察0~Nyquist频率范围的幅频特性就可以了。
3、利用FFT进行功率谱的噪声分析
已知带有测量噪声信号其中f1=50Hz，f2=120Hz，为均值为零、方差为1的随机信号，采样频率为1000Hz，数据点数N=512。试绘制信号的功率谱图.
MATLAB程序如下
t=0：0。001：0。6；
x=sin（2＊pi＊50*t）+sin(2*pi*120*t);
y=x+2＊randn（1,length（t))；
Y=fft(y,512）;
P=Y。＊conj（Y）/512; %求功率
f=1000*（0:255)/512;
subplot(2，1,1）；
plot(y)；
subplot（2，1，2）;
plot(f,P（1：256))；
程序运行结果如图6—2

图6-2
4、序列长度和FFT的长度对信号频谱的影响。
已知信号
其中f1=15Hz，f2=40Hz，采样频率为100Hz。
在下列情况下绘制其幅频谱。
Ndata=32，Nfft=32；
Ndata=32,Nfft=128;
MATLAB程序如下
fs=100；
Ndata=32; Nfft=32;
n=0：Ndata—1；
t=n/fs;
x=0。5＊sin(2*pi*15＊t)+2＊sin（2*pi*40*t）；
y=fft（x，Nfft）；
mag=abs(y)；
f=(0:length（y)—1）'＊fs/length(y）;
subplot(2,1,1）
plot(f（1:Nfft/2),mag(1:Nfft/2))
title('Ndata=32,Nfft=32')
程序运行结果如图6-3所示

图6-3
5、快速卷积的FFT算法
在MATLAB实现卷积的函数为CONV，对于N值较小的向量，这是十分有效的。对于N值较大的向量卷积可用FFT加快计算速度.
由DFT性质可知,若DFT[x1（n)]=X1（k),DFT[x2(n）］=X2(k)则

若DFT和IDFT均采用FFT和IFFT算法,可提高卷积速度。
五、实验任务
1、输入并运行例题程序，熟悉基本指令的使用。
2、比较定义式计算傅里叶变换和用快速算法计算傅里叶变换所用时间。
3、比较卷积函数与快速卷积运算所用时间。
（提示：clock函数读取瞬时时钟
etime(t1,t2)函数计算时刻t1，t2间所经历的时间。）
六、实验报告
1、简述实验目的、原理.
2、写出上机调试通过的实验任务的程序并比较它们运行时间的优劣.
3、思考题：通过实验总结快速算法的优越性。

实验七 模拟原型滤波器设计
一、实验目的
1、加深对模拟滤波器基本类型、特点及其主要性能指标的掌握。
2、掌握模拟低通滤波器原型的设计方法。
3、掌握MATLAB工具箱函数的调用.
二、实验原理
1、模拟滤波器
输入信号和输出信号均为连续信号，冲激响应也是连续的滤波器,称为模拟滤波器。
模拟滤波器从功能上可以分为低通、高通、带通、带阻以及全通滤波器。
实际使用中理想滤波器是不可实现的,必须设计一个因果可实现的滤波器去逼近。通常,通带和阻带都允许存在一定的误差范围，即通带不一定是完全水平的，阻带也不一定绝对衰减到零。在通带和阻带之间允许设置一定宽度的过渡带。
2、典型的模拟滤波器
典型的模拟滤波器有巴特沃斯滤波器、切比雪夫滤波器(一型/二型）、椭圆滤波器等。每种滤波器都有其不同的特点。
巴特沃斯滤波器具有单调下降的幅频特性，通带和阻带幅频都比较平坦.
切比雪夫1型滤波器在通带内具有等波动的幅频特性。
切比雪夫2型滤波器在阻带内具有等波动的幅频特性.
椭圆滤波器在通带和阻带内均具有等波动的幅频特性.
三、实验用函数
1、buttord
功能：确定巴特沃斯滤波器的最小阶数和截止频率。
调用格式：
［n,wn]=buttord(wp，ws，rp，rs)：计算巴特沃斯数字滤波器的阶数和截止频率。
[n,wn］=buttord（wp，ws，rp，rs，’s’):计算巴特沃斯模拟滤波器的阶数和截止频率。
其中：wp通带截止频率，ws阻带截止频率，rp通带衰减，rs阻带衰减.wp，ws为一元向量时，为低通或高通滤波器,wp，ws为二元向量时，为带通或带阻滤波器。
2、cheb1ord
功能：确定切比雪夫1型滤波器的最小阶数和通带截止频率。
调用格式:
［n,wn］=cheb1ord(wp,ws，rp，rs）：计算切比雪夫1型数字滤波器的最小阶数和通带截止频率.
[n，wn]= cheb1ord (wp，ws,rp,rs，’s’）：计算切比雪夫1型模拟滤波器的最小阶数和通带截止频率。
3、cheb2ord
功能：确定切比雪夫2型滤波器的最小阶数和阻带截止频率.
调用格式：
［n，wn]=cheb2ord(wp,ws，rp，rs):计算切比雪夫2型数字滤波器的最小阶数和阻带截止频率。
[n,wn]= cheb2ord （wp,ws,rp，rs，’s’）：计算切比雪夫2型模拟滤波器的最小阶数和阻带截止频率。
4、ellipord
功能：确定椭圆滤波器的最小阶数和通带截止频率。
调用格式：
[n，wn］=ellipord(wp,ws，rp,rs)：计算椭圆数字滤波器的最小阶数和通带截止频率.
［n，wn]= ellipord(wp,ws，rp，rs，’s’)：计算椭圆模拟滤波器的最小阶数和通带截止频率。
5、buttap
功能：巴特沃斯模拟低通滤波器原型。
调用格式：
［z,p，k]=buttap(n)：设计巴特沃斯模拟低通滤波器原型，其传递函数为

此时z为空阵。巴特沃斯滤波器由通带内最平坦、总体上单调的幅度特性来表征。
6、cheb1ap
功能：切比雪夫1型模拟低通滤波器原型。
调用格式：
［z,p，k］=cheb1ap（n,rp）:设计切比雪夫1型模拟低通滤波器原型，其通带内的波纹系数为rp分贝,传递函数为

此时z为空阵。切比雪夫1型滤波器为通带内等波纹、阻带内单调的滤波器,其极点均匀分布在左半平面的椭圆上.
7、cheb2ap
功能：切比雪夫2型模拟低通滤波器原型.
调用格式：
［z，p,k]=cheb2ap（n，rs）：设计切比雪夫2型模拟低通滤波器原型，其阻带内的波纹系数小于rs分贝，传递函数为

切比雪夫2型滤波器为通带内单调、阻带内等波纹的滤波器，其极点位置为cheb1ap极点位置的倒数.
8、ellipap
功能：椭圆模拟低通滤波器原型。
调用格式：
［z，p,k]=ellipap(n，rp,rs)：设计椭圆模拟低通滤波器原型，其通带内的波纹系数为rp分贝，阻带内的波纹系数小于通带的rs分贝,传递函数为

椭圆滤波器为通带内和阻带内等波纹的滤波器，它具有比巴特沃斯和切比雪夫更陡的下降斜率,但会损失通带和阻带的波纹指标。
四、实例
通过模拟滤波器原型设计一个巴特沃斯模拟低通滤波器，要求通带截止频率fp=2kHz，通带最大衰减Rp<=1dB，阻带截止频率fs=5kHz，阻带最小衰减As〉=20dB.
MATLAB程序:
fp=2000；
fs=5000;
rp=1;as=20；
［n，wn]=buttord（fp,fs,rp,as,’s’）;
[z,p,k］=buttap(n)；
[b,a]=zp2tf(z,p,k)；
freqs(b,a)
程序运行结果如图7-1所示。

图7—1
五、实验内容
1、输入并运行例题程序，熟悉基本指令的使用。
2、设计一个模拟原型低通滤波器，要求通带截止频率fp=6kHz,通带最大衰减Rp〈=1dB，阻带截止频率fs=15kHz，阻带最小衰减As>=30dB.要求：分别利用巴特沃斯、切比雪夫、椭圆等滤波器来实现.熟悉几种经典滤波器的基本使用。
六、实验报告
1、简述实验目的、原理。
2、写出上机调试通过的实验任务的程序并描述其图形曲线.

实验八 数字滤波器设计-—IIR
一、实验目的
1、加深对数字滤波器基本类型、特点及其主要性能指标的掌握。
2、掌握IIR数字滤波器的设计方法。
3、掌握MATLAB工具箱函数的调用。
二、实验原理
1、数字滤波器
数字滤波是数字信号处理技术的重要内容。和模拟滤波器一样，数字滤波器的主要功能是对数字信号进行处理，保留数字信号中的有用成分,去除信号中的无用成分.
数字滤波器是具有一定传输特性的数字信号处理装置。它的输入和输出均为离散的数字信号，借助数字器件或一定的数值计算方法，对输入信号进行