数字信号处理实验指导书

数字信号处理 实 验 指 导 书

0 / 49

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

一、实验目的

1、了解常用时域离散信号及其特点； 2、掌握MATLAB程序的编程方法； 3、熟悉MATLAB函数的调用方法。

二、实验原理

在时间轴上的离散点取值的信号，称为离散时间信号。离散时间信号只在某些离散的瞬时给出函数值，而在其他时刻无定义。它是时间上不连续按一定先后次序排列的一组数的集合，称为时间序列，用x(n)表示，n取整数代表时间的离散时刻。

在MATLAB中用向量来表示一个有限长度的序列。 常用离散信号： 1、单位抽样序列

(n)1n=n01n0 或(n-n0)=0n00nn02、单位阶跃序列

1nn01n0u(n)或u(n-n0)

0n00nn03、实指数序列

x(n)an

4、复指数序列

x(n)e(j)t

5、正(余)弦序列

x(n)Umsin(0n)

6、随机序列

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

7、周期序列

x(n)x(nN)1 / 49

三、实验用函数

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);

2 / 49

stem(n,x)

程序运行结果如图1-1所示：

图1-1 单位抽样序列

例1.2实数指数序列（运算符“.^”）

Matlab程序如下： n=[0:10]; x=0.9.^n; stem(n,x)

程序运行结果如图1-2所示

图1-2 实数指数序列

3 / 49

例1.3复数指数序列（x(n)eMatlab程序如下：

(0.1j0.3)n(10n10)）

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正、余弦序列（x(n)Umsin(0n)）

Matlab程序如下： n=[0:10];

x=3*cos(0.1*pi*n+pi/3); stem(n,x)

程序运行结果如图1-4所示

4 / 49

图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所示

5 / 49

图1-5 周期序列

五、实验任务

1、调试部分例题程序，掌握Matlab基本操作方法。 2、编写程序，完成下列函数波形： 1）利用zeros函数生成单位抽样序列；

2）利用zeros函数和ones函数生成单位阶跃序列；

六、实验报告

1、简述实验目的、原理。

2、写出上机调试通过的实验任务的程序并描述其图形曲线。

6 / 49

实验二 离散序列的基本运算

一、实验目的

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;

7 / 49

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）序列相乘

序列相加是指两个序列中相同序号的序列值逐项对应相乘，形成新的序列。 参加运算的两个序列的维数不同时处理方法与序列相加相同。

8 / 49

五、实验任务

1、理解序列运算的性质，了解函数语句的意义。 2、利用例题函数完成下列序列运算 1）已知x1(n)=u(n+1) (-3x2(n)=u(n-3) (-4求：x(n)=x1(n)+x2(n)

2）已知x1(n)=3e-0.25n (-2六、实验报告

1、简述实验目的和原理。

2、列写上机调试通过的程序，并描绘其波形曲线。

9 / 49

实验三 离散卷积的原理及应用

一、实验目的

1、通过实验进一步理解卷积定理，了解卷积过程； 2、掌握应用线性卷积求解离散时间系统响应的基本方法。

二、实验原理

对于线性移不变离散系统，任意的输入信号x(n)可以用(n)及其位移的线性组合来表示，即

x(n)kx(k)(nk)

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

y(n)kx(k)h(nk)

称为离散系统的线性卷积，简写为

y(n)=x(n)*h(n)

也就是说，如果已知系统的冲激响应，将输入信号与系统的冲激响应进行卷积运算，即可求得系统的响应。

三、实验用函数

1、卷积函数conv

功能：进行两个序列的卷积运算。 调用格式：

y=conv(x,h)；用于求解两有限长序列的卷积。 2、sum

功能：求各元素之和。 调用格式：

y=sum(x)；求序列x中各元素之和。 3、hold

功能：控制当前图形窗口是否刷新的双向切换开关。 调用格式：

hold on：使当前图形窗口中的图形保持且不被刷新，准备接受绘制新的图形。 hold off：使当前图形窗口中的图形不具备不被刷新的性质。

10 / 49

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数字低通滤波器的系统函数公式为

0.13210.3963z10.3963z20.1321z3H(z)

10.34319z10.60439z20.20407z3输入一个矩形信号序列 x=square(n/5) (-211 / 49

六、实验报告

1、简述实验目的、原理。

2、写出上机调试通过的实验任务的程序并描述其图形曲线。

12 / 49

实验四 离散傅立叶级数

一、实验目的

1、加深对离散周期序列傅里叶级数基本概念的理解。

2、掌握MATLAB求解周期序列傅里叶级数变换和逆变换的方法。 3、观察离散周期序列的重复周期数对频谱特性的影响。

二、实验原理

离散时间序列x(n)满足x(n)=x(n+rN)，称为离散周期序列，其中N为周期，x(n)为主值序列。

周期序列可用离散傅里叶级数表示成

1x(n)NX(k)ek0N1j2knNIDFS[X(k)] n=0,1,…,N-1

其中，X(k)是周期序列离散傅里叶级数第K次谐波分量的系数，也称为周期序列的频谱，可表示为

X(k)x(n)en0N1j2knNDFS[x(n)] k=0,1,…,N-1

上面两式是周期序列的一对傅里叶级数变换对。

令WNej2N，以上两式可简写为：

nkX(k)DFS[x(n)]x(n)WNn0N11x(n)IDFS[X(k)]NX(k)Wn0N1

nkN三、实验用函数

1、mod

功能：模除求余。 调用格式：

mod(x,m)：x整除m取正余数。 2、floor

功能：向-舍入为整数。

13 / 49

调用格式：

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点。用傅里叶级数求信号的幅度和相位频谱；求傅里叶级数逆变换的图形，与原信号图形进行比较。

14 / 49

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无限增大，研究其极限情况。这样可以观察信号序列由非周期到周期

15 / 49

变换时，频谱由连续谱逐渐向离散谱过渡的过程。

16 / 49

例：已知一个矩形序列的脉冲宽度占整个周期的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

16 / 49

图4-2

从上图可以看出，信号序列的周期数越多，则频谱越是向几个频点集中，当信号周期数趋于无穷大时 ，频谱转化为离散谱。

五、实验任务

1、输入并运行例题程序，熟悉基本指令的使用。

2、已知一个信号序列的主值为x(n)=[0,1,2,3,2,1,0]，显示两个周期的信号序列波形。 求解：用傅里叶级数求信号的幅度和相位频谱；求傅里叶级数逆变换的图形，与原信号图形进行比较。

六、实验报告

1、简述实验目的、原理。

2、写出上机调试通过的实验任务的程序并描述其图形曲线。

17 / 49

实验五 离散傅里叶变换

一、实验目的

1、加深对离散傅里叶变换基本概念的理解。

2、了解有限长序列傅里叶变换与周期序列傅里叶级数的联系。 3、熟悉相关函数的使用方法。

二、实验原理

有限长序列的傅里叶变换和逆变换

对于非周期序列，在实际中常常使用有限长序列。有限长序列x(n)表示为

x(n)0nN1x(n)

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

X(k)DFT[x(n)]x(n)en0N1j2nkN

1x(n)IDFT[x(k)]NDFT与DFS的关系

X(k)ek0N1j2knN

比较两者的变换对，可以看出两者的区别仅仅是将周期序列换成了有限长序列。 有限长序列x(n)可以看作是周期序列x(n)的一个周期；反之周期序列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)]图形。

18 / 49

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-119 / 49

由上图可看出，与周期序列不同的是，有限长序列本身是仅有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：

20 / 49

图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、写出上机调试通过的实验任务的程序并描述其图形曲线。

21 / 49

实验六 快速傅里叶变换

一、实验目的

1、加深对快速傅里叶变换基本理论的理解。 2、了解用MATLAB语言进行快速傅里叶变换的方法。 3、掌握常用函数的调用方法。

二、实验原理

DFT是在时域和频域均为离散序列的变换方法，它适用于有限长序列。但如果按照变换公式进行运算的话，当序列长度很大时，将占用很大的内存空间，运算时间也会很长，无法实时处理问题。

快速傅里叶变换是用于提高DFT运算的高速运算方法的统称，FFT是其中的一种，FFT不是一种新的变换形式，它仅仅只是一种快速算法。FFT主要有时域抽取算法和频域抽取算法，基本思想是将一个长度为N的序列分解成多个短序列再进行运算，如基2算法、基4算法等等，从而可以大大缩短运算时间。

三、实验用函数

1、fft

功能：一维基2快速傅里叶变换。 调用格式：

y=fft(x)：利用FFT算法计算矢量x的离散傅里叶变换。当x为2的幂次方时，采用高

22 / 49

速基2FFT算法，否则为稍慢的混合算法。

23 / 49

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

23 / 49

-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

24 / 49

图6-1

如图所示，由于信号采样频率为100Hz，故其莱奎斯特频率为50Hz，图中整个频谱图是以莱奎斯特频率为轴对称的。因此利用FFT对信号作谱分析时，只要考察0~Nyquist频率范围的幅频特性就可以了。

3、利用FFT进行功率谱的噪声分析

已知带有测量噪声信号x(t)sin(2f1t)sin(2f2t)2(t)其中f1=50Hz，f2=120Hz，(t)为均值为零、方差为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);

25 / 49

plot(f,P(1:256)); 程序运行结果如图6-2

图6-2

4、序列长度和FFT的长度对信号频谱的影响。 已知信号x(t)0.5sin(2f1t)2sin(2f2t) 其中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);

26 / 49

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)则

x1(n)*x2(n)IDFT[X1(k)•X2(k)]IDFT[DFT[x1(n)]•DFT[x2(n)]] 若DFT和IDFT均采用FFT和IFFT算法，可提高卷积速度。

五、实验任务

1、输入并运行例题程序，熟悉基本指令的使用。

2、比较定义式计算傅里叶变换和用快速算法计算傅里叶变换所用时间。 3、比较卷积函数与快速卷积运算所用时间。 （提示：clock函数读取瞬时时钟

etime(t1,t2)函数计算时刻t1,t2间所经历的时间。）

六、实验报告

1、简述实验目的、原理。

27 / 49

2、写出上机调试通过的实验任务的程序并比较它们运行时间的优劣。 3、思考题：通过实验总结快速算法的优越性。

实验七 模拟原型滤波器设计

一、实验目的

1、加深对模拟滤波器基本类型、特点及其主要性能指标的掌握。 2、掌握模拟低通滤波器原型的设计方法。 3、掌握MATLAB工具箱函数的调用。

二、实验原理

1、模拟滤波器

输入信号和输出信号均为连续信号，冲激响应也是连续的滤波器，称为模拟滤波器。 模拟滤波器从功能上可以分为低通、高通、带通、带阻以及全通滤波器。

实际使用中理想滤波器是不可实现的，必须设计一个因果可实现的滤波器去逼近。通常，通带和阻带都允许存在一定的误差范围，即通带不一定是完全水平的，阻带也不一定绝对衰减到零。在通带和阻带之间允许设置一定宽度的过渡带。

2、典型的模拟滤波器

28 / 49

典型的模拟滤波器有巴特沃斯滤波器、切比雪夫滤波器（一型/二型）、椭圆滤波器等。每种滤波器都有其不同的特点。

巴特沃斯滤波器具有单调下降的幅频特性，通带和阻带幅频都比较平坦。 切比雪夫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

29 / 49

功能：确定椭圆滤波器的最小阶数和通带截止频率。 调用格式：

[n,wn]=ellipord(wp,ws,rp,rs)：计算椭圆数字滤波器的最小阶数和通带截止频率。 [n,wn]= ellipord(wp,ws,rp,rs,’s’)：计算椭圆模拟滤波器的最小阶数和通带截止频率。

5、buttap

功能：巴特沃斯模拟低通滤波器原型。 调用格式：

[z,p,k]=buttap(n)：设计巴特沃斯模拟低通滤波器原型，其传递函数为

Ha(s)k

(sp(1))(sp(2))...(sp(n))此时z为空阵。巴特沃斯滤波器由通带内最平坦、总体上单调的幅度特性来表征。 6、cheb1ap

功能：切比雪夫1型模拟低通滤波器原型。 调用格式：

[z,p,k]=cheb1ap(n,rp)：设计切比雪夫1型模拟低通滤波器原型，其通带内的波纹系数为rp分贝，传递函数为

Ha(s)k

(sp(1))(sp(2))...(sp(n))此时z为空阵。切比雪夫1型滤波器为通带内等波纹、阻带内单调的滤波器，其极点均匀分布在左半平面的椭圆上。

7、cheb2ap

功能：切比雪夫2型模拟低通滤波器原型。 调用格式：

[z,p,k]=cheb2ap(n,rs)：设计切比雪夫2型模拟低通滤波器原型，其阻带内的波纹系数小于rs分贝，传递函数为

Ha(s)k(sz(1))(sz(2))...(sz(n))

(sp(1))(sp(2))...(sp(n))切比雪夫2型滤波器为通带内单调、阻带内等波纹的滤波器，其极点位置为cheb1ap

30 / 49

极点位置的倒数。

31 / 49

8、ellipap

功能：椭圆模拟低通滤波器原型。 调用格式：

[z,p,k]=ellipap(n,rp,rs)：设计椭圆模拟低通滤波器原型，其通带内的波纹系数为rp分贝，阻带内的波纹系数小于通带的rs分贝，传递函数为

Ha(s)k(sz(1))(sz(2))...(sz(n))

(sp(1))(sp(2))...(sp(n))椭圆滤波器为通带内和阻带内等波纹的滤波器，它具有比巴特沃斯和切比雪夫更陡的下降斜率，但会损失通带和阻带的波纹指标。

四、实例

通过模拟滤波器原型设计一个巴特沃斯模拟低通滤波器，要求通带截止频率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所示。

31 / 49

图7-1

五、实验内容

1、输入并运行例题程序，熟悉基本指令的使用。

2、设计一个模拟原型低通滤波器，要求通带截止频率fp=6kHz，通带最大衰减Rp<=1dB，阻带截止频率fs=15kHz，阻带最小衰减As>=30dB。要求：分别利用巴特沃斯、切比雪夫、椭圆等滤波器来实现。熟悉几种经典滤波器的基本使用。

六、实验报告

1、简述实验目的、原理。

2、写出上机调试通过的实验任务的程序并描述其图形曲线。

实验八 数字滤波器设计——IIR

一、实验目的

1、加深对数字滤波器基本类型、特点及其主要性能指标的掌握。 2、掌握IIR数字滤波器的设计方法。 3、掌握MATLAB工具箱函数的调用。

二、实验原理

32 / 49

1、数字滤波器

数字滤波是数字信号处理技术的重要内容。和模拟滤波器一样，数字滤波器的主要功能是对数字信号进行处理，保留数字信号中的有用成分，去除信号中的无用成分。

数字滤波器是具有一定传输特性的数字信号处理装置。它的输入和输出均为离散的数字信号，借助数字器件或一定的数值计算方法，对输入信号进行处理，改变输入信号的波形或频谱，达到保留信号中有用成分去除无用成分的目的。如果加上A/D、D/A转换，则可以用于处理模拟信号。

2、数字滤波器分类

按时域特性，数字滤波器可以分为无限冲激响应数字滤波器（IIR）和有限长冲激响应数字滤波器（FIR）两类。

按频域特性，数字滤波器和模拟滤波器一样可分为低通、高通、带通、带阻等。数字滤波器是一个离散时间系统，其频率特性具有周期性，因此我们在讨论的频率范围仅在

0范围内，相应的标准化频率在0~1之间。

3、IIR数字滤波器的设计方法

由于它的脉冲相应序列是无限的，故称为无限冲激响应滤波器。IIR滤波器的设计是根据滤波器的性能指标要求，设计滤波器的分子和分母多项式。IIR滤波器常用设计方法如下表所示：

设计方法 说明 首先设计满足性能要求的模拟滤波器，再离散化为数字滤波器 完全设计函数 butter、cheby1、cheby2、ellip 阶数估计函数 buttord、cheb1ord、cheb2ord、ellipord 经典设计法 低通模拟原型滤波器函数 buttap、cheb1ap、cheb2ap、ellipord 频率转换函数 lp2lp、lp2hp、lp2bp、lp2sp 滤波器离散化函数 工具函数 33 / 49

在离散域内用最小二乘法逼直接设计法 近给定的幅频特性 采用参数模型逼近给定时域或频率响应求得滤波器 参数模型设计法 bilinear、impinvar yulewalk 时域建模函数 lpc、prony、stmcb 频域建模函数 invfreqs(模拟)、invfreqz(数字) 设计一般化低通模拟滤波器，maxflat 最大平滑滤波器 其零点多余极点 （1）IIR滤波器经典设计法 1、模拟滤波器变换法

首先根据滤波器的技术指标设计相应的模拟滤波器，然后再将设计好的模拟滤波器变换成满足给定指标的数字滤波器。

具体设计步骤为：

1）对设计性能指标中的频率指标进行转换使其满足模拟滤波器原型设计性能指标； 2）估计模拟滤波器最小阶数和边界频率；Matlab提供的函数（buttord,cheb1ord, cheb2ord,ellipord）；

3）设计模拟低通滤波器原型； Matlab提供的函数（buttap,cheb1ap,cheb2ap,ellipap）；

4）由模拟低通原型经频率变换获得模拟滤波器； Matlab提供的函数（lp2lp,lp2hp,lp2bp, lp2bs）

5）将模拟滤波器离散化获得IIR数字滤波器。 Matlab提供的函数（bilinear,impinvar） 2、IIR滤波器完全设计函数

MATLB工具箱还提供了IIR滤波器设计的完全工具函数，用户只要调用这些工具函数即可一次性完成设计，而不需要分步进行实现。在使用这些工具函数时需要注意其频率要采用标准化频率。

34 / 49

（2）IIR滤波器直接设计法

经典设计法只限于几种标准的低通、高通、带通、带阻滤波器，而对于具有任意形状或多频带滤波器的设计是无能为力的。

采用最小二乘法拟合给定的幅频响应，使设计的滤波器幅频特性逼近期望的频率特性，这种方法称为IIR滤波器的直接设计法。

三、实验用函数

1、频率变换函数（lp2lp、lp2hp、lp2bp、lp2bs）

功能：模拟低通原型到低通、高通、带通、带阻模拟滤波器变换。 调用格式：

[bt,at]= lp2lp/lp2hp/lp2bp/lp2bs(b,a,wn)（wn为滤波器截至频率） 2、impinvar

功能：用脉冲响应不变法实现模拟到数字的滤波器变换。 调用格式：

[bd,ad]=impinvar(b,a,fs)：将模拟滤波器系数b、a变换为数字滤波器系数bd、ad，两者的冲激响应不变。

3、bilinear

功能：将s域映射到z域的标准方法。 调用格式：

[numd,dend]=bilinear(num,den,fs)：将模拟域传递函数变换为数字域传递函数。 [numd,dend]=bilinear(num,den,fs,fp)：将模拟域传递函数变换为数字域传递函数，fs为取样频率，fp为通带截止频率。

[zd,pd,kd]=bilinear(z,p,k,fs)：将模拟域零极点增益系数变换到数字域，fs为取样频率。

[zd,pd,kd]=bilinear(z,p,k,fs,fp)：将模拟域零极点增益系数变换到数字域，fs为取样频率，fp为通带截止频率。

[ad,bd,cd,dd]=bilinear(a,b,c,d,fs)：将模拟域状态变量系数变换到数字域，fs为取样频率。

4、butter

功能：巴特沃斯模拟和数字滤波器设计。

35 / 49

调用格式：

[b,a]=butter(n,wn,’ftype’,’s’) ftype为滤波器类型：

‘high’为高通滤波器，截止频率wn

‘stop’为带阻滤波器，截止频率wn=[w1,w2] 缺省时为低通或带通滤波器 5、cheby1

功能：切比雪夫1型模拟和数字滤波器设计。 调用格式：

[b,a]=cheby1(n,rp,wn,’ftype’,’s’)。 6、cheby2

功能：切比雪夫2型模拟和数字滤波器设计。。 调用格式：

[b,a]=cheby2(n,rs,wn,’ftype’,’s’)。 7、ellip

功能：椭圆模拟和数字滤波器设计。 调用格式：

[b,a]=ellip(n,rp,rs,wn,’ftype’,’s’)。

四、实例

例1、用脉冲响应不变法设计Butterworth低通数字滤波器，要求通带频率为

00.2，通带波纹小于1dB，阻带在内0.3，幅度衰减大于15dB，采样周

期为Ts=0.01s。

MATLAB程序

wp=0.2*pi; ws=0.3*pi; rp=1;rs=15;ts=0.01;Nn=128; Wp=wp/ts; Ws=ws/ts;

[N,Wn]=buttord(Wp,Ws,rp,rs,'s'); [z,p,k]=buttap(N); [Bap,Aap]=zp2tf(z,p,k);

36 / 49

[b,a]=lp2lp(Bap,Aap,Wn); [bz,az]=impinvar(b,a,1/ts); freqz(bz,az,Nn,1/ts) 程序运行结果如图8-1所示

图8-1

例2、设计一个Butterworth高通数字滤波器，通带边界频率为300Hz，阻带边界频率为200Hz，通带波纹小于1dB，阻带衰减大于20dB，采样频率为1000Hz。

MATLAB程序 fs=1000;

wp=300/(fs/2); ws=200/(fs/2); rp=1; rs=15; Nn=128; [N,Wn]=buttord(wp,ws,rp,rs); [b,a]=butter(N,Wn,'high'); freqz(b,a,Nn,fs) 程序运行结果如图8-2所示

37 / 49

图8-2

五、实验内容

1、输入并运行例题程序，熟悉基本指令的使用。

2、用脉冲响应不变法设计一个切比雪夫型数字带通滤波器，要求通带频率

0.30.7，通带最大衰减Rp=1dB，阻带截止频率s10.1,s20.9，阻带最小

衰减As=15dB，滤波器采样频率Fs=2000Hz。

六、实验报告

1、简述实验目的、原理。

2、写出上机调试通过的实验任务的程序并描述其图形曲线。

38 / 49

实验九 数字滤波器设计——FIR

一、实验目的

1、加深对数字滤波器基本类型、特点及其主要性能指标的掌握。 2、掌握FIR数字滤波器的设计方法。 3、掌握MATLAB工具箱函数的调用。

二、实验原理

1、数字滤波器

数字滤波是数字信号处理技术的重要内容。和模拟滤波器一样，数字滤波器的主要功能是对数字信号进行处理，保留数字信号中的有用成分，去除信号中的无用成分。

数字滤波器是具有一定传输特性的数字信号处理装置。它的输入和输出均为离散的数字信号，借助数字器件或一定的数值计算方法，对输入信号进行处理，改变输入信号的波形或频谱，达到保留信号中有用成分去除无用成分的目的。如果加上A/D、D/A转换，则可以用于处理模拟信号。

2、数字滤波器分类

按时域特性，数字滤波器可以分为无限冲激响应数字滤波器（IIR）和有限长冲激响应数字滤波器（FIR）两类。

按频域特性，数字滤波器和模拟滤波器一样可分为低通、高通、带通、带阻等。数字滤波器是一个离散时间系统，其频率特性具有周期性，因此我们在讨论的频率范围仅在

0范围内，相应的标准化频率在0~1之间。

3、FIR数字滤波器的设计方法

由于它的脉冲相应序列是有限长，故称为有限长冲激响应滤波器。FIR滤波器与IIR滤波器相比较其突出优点是，在保证满足滤波器幅频响应要求的同时，还可以获得严格的线性相位特性，这对于高保真的信号处理。如语音处理、数据处理和测试等是十分重要的。它的主要缺点是，达到相同性能指标所需滤波器阶数要高得多，延迟也比较大。FIR滤波器常用设计方法如下表所示：

设计方法 窗函数设计法 kaiserord 说明 理想滤波器加窗处理 fir1、fir2 工具函数 39 / 49

最小平方误差最小化逼近理最优化设计法 想幅频响应 约束的最小二乘逼近 任意响应设计 包括复响应和非线性相位 具有光滑、正弦过渡带的低通升余弦函数 滤波器设计 （1）FIR滤波器窗函数设计法

采用参数模型逼近给定时域或频率响应求得滤波器 firls、remez、remezord fircls fircls1 具有任意响应的FIR滤波器，cremez firrcos FIR滤波器设计的主要任务是根据给定的性能指标确定滤波器的系数b。即系统单位脉冲序列h(n)，它是一个有限长序列。

具体设计步骤为：

1）对滤波器理想特性Hd(e)求出其单位脉冲响应hd(n)；

对于理想的数字低通滤波器频率响应，可以通过下面子程序来实现： function hd=ideal_lp(w,N) %hd为理想脉冲响应序列 %w为截至频率 %N为理想滤波器的长度 alpha=(N-1)/2; n=[0:(N-1)];

m=n-alpha+eps; %加一个极小值以避免0作为除数 hd=sin(w*m)./(pi*m);

2）由性能指标确定窗函数W(n)和窗口长度N，由过渡带宽度近似于窗函数主瓣宽求得窗口长度N。

3）求得实际滤波器的单位脉冲序列h(n)=hd(n)*w(n) h(n)即为所设计FIR滤波器系数向量。 （2）最优FIR滤波器设计法

jw

40 / 49

MATLAB工具箱提供了比基于窗函数法的FIR滤波器设计工具函数更通用的函数firls、remez。它们采用不同的优化方法设计最优的标准的和多频带FIR数字滤波器。

函数firls是函数fir1和fir2的扩展，其基本设计准则是利用最小二乘法使期望的幅频响应和实际的频率响应之间的整体误差最小。

函数remez实现Parks-McClellan算法，这种算法利用Remez交换算法和Chebyshev近似理论来设计滤波器，使实际频率响应拟合期望频率响应达到最优。从实际和理想频响之间最大误差最小化的观点来看，函数remez设计的滤波器是最优的，因此，有时称之为最小滤波器。在频域内，滤波器呈现等波纹特点，因此也称为等波纹滤波器。该设计方法是FIR滤波器设计中最流行的和应用最广的。

三、实验用函数

1、boxcar 功能：矩形窗。 调用格式：

w=boxcar(n)：产生一个长度为n的矩形窗函数。 2、triang 功能：三角窗。 调用格式：

w=triang (n)：得到n点的三角窗函数。 3、bartlett 功能：巴特利特窗。 调用格式：

w=bartlett(n)：得到n点的巴特利特窗函数。 4、hamming 功能：哈明窗。 调用格式：

w=hamming(n)：得到n点的哈明窗函数。 5、hanning 功能：汉宁窗。 调用格式：

41 / 49

w=hanning(n)：得到n点的汉宁窗函数。 6、blackman 功能：布莱克曼窗。 调用格式：

w=blackman(n)：得到n点的布莱克曼窗函数。 7、chebwin

功能：切比雪夫窗。 调用格式：

w=chebwin(n)：得到n点的切比雪夫窗函数。 8、kaiser 功能：凯塞窗。 调用格式：

w=kaiser(n,beta)：得到n点的凯塞窗函数，其中beta为影响窗函数旁瓣的参数。 9、fir1

功能：基于窗函数的FIR数字滤波器设计——标准频率响应，以经典方法实现加窗线性相位FIR滤波器设计，可设计出标准的低通、高通、带通、带阻滤波器。

调用格式：

b=fir1(n,wn,’ftype’,window)：设定不同参数设计各种加窗的滤波器。 10、fir2

功能：设计具有任意形状频率响应的FIR数字滤波器。 调用格式：

b=fir2(n,f,m,npt,lap,window)：设定不同参数设计各种加窗的滤波器。

其中：f,m分别为滤波器期望幅频响应的频率向量和幅值向量，取值在0~1之间；npt为对频率响应进行内插点数，缺省为512,；lap定义一个区域尺寸，函数fir2在重复频率点周围建立这个区域并提供光滑、陡峭的过渡频率响应，缺省值为25；b为FIR滤波器系数向量。

11、firls、remez 功能：最优化FIR滤波器。

42 / 49

调用格式：

b=firls(n,f,a)、b=remez(n,f,a)：基本形式最优滤波器，f,a分别是理想频率特性频率向量和幅值向量。

b=firls(n,f,a,w)、b=remez(n,f,a,w)：加权最优滤波器，w为权向量，为f和a向量长度的一半，一个频带必须对应一个权值。通过设置权值使不同频段的误差最小化得到不同程度重视。

四、实例

例1、用矩形窗设计一个FIR数字低通滤波器，要求：N=64，截止频率为c0.4，描绘理想和实际滤波器的脉冲响应、窗函数及滤波器的幅频响应曲线。

MATLAB程序

wc=0.4*pi; N=64; n=0:N-1;

hd=ideal_lp(wc,N); windows=(boxcar(N))'; b=hd.*windows;

subplot(3,1,1); stem(n,hd); title('理想脉冲响应'); subplot(3,1,2); stem(n,windows); title('窗函数特性'); subplot(3,1,3); stem(n,b); title('实际脉冲响应'); 程序运行结果如图9-1所示

43 / 49

图9-1

例2、用MATLAB信号处理箱提供的fir1函数，设计一个FIR数字高通滤波器。要求：通带截止频率为fp=450Hz，Rp=0.5dB；阻带截止频率为fs=300Hz，As=20dB；采样频率Fs=2000Hz。描绘滤波器的脉冲响应、窗函数及滤波器的幅频响应和相频响应曲线。

求解频率特性的工具函数

function [db,mag,pha,grd,w] = freqz_m(b,a); [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); grd = grpdelay(b,a,w); MATLAB程序

fs=300; fp=450; Fs=2000; ws=fs/(Fs/2)*pi; wp=fp/(Fs/2)*pi; %归一化

deltaw=wp-ws; %计算过渡带的宽度

N0=ceil(6.1*pi/deltaw); %计算滤波器窗口长度

N=N0+mod(N0+1,2); %为实现FIR类型I偶对称滤波器，应确保N为奇数 n=0:N-1;

windows=triang(N);

wc=(ws+wp)/2/pi; %截止频率取归一化通阻带频率的平均值 b=fir1(N-1,wc,'high',windows); %求系统函数系数 [db,mag,pha,grd,w]=freqz_m(b,1); %求解频率特性 subplot(2,2,1); stem(n,b); subplot(2,2,2); stem(n,windows); subplot(2,2,3); plot(w/pi,db); subplot(2,2,4); plot(w/pi,pha);

44 / 49

程序运行结果如图9-2所示

图9-2

五、实验内容

1、输入并运行例题程序，熟悉基本指令的使用。 2、用窗函数法设计FIR数字带通滤波器

选择合适的窗函数用工具箱提供的fir1函数设计一个FIR数字带通滤波器，要求：下阻带截止频率fs1100Hz,As65dB，通带低端截止频率fp1150Hz,Rp0.05dB；通带高端端截止频fp2350Hz,Rp0.05dB；上阻带截止频率fs2400Hz,As65dB；描绘滤波器的脉冲响应、窗函数及滤波器的幅频响应和相频响应曲线。。

六、实验报告

1、简述实验目的、原理。

2、写出上机调试通过的实验任务的程序并描述其图形曲线。

注：合同范本有风险，使用需谨慎，法律是经验性极强的领域，范本无法思考和涵盖全面，最好找专业律师起草或审核后使用，谢谢您的关注！

45 / 49