西科大《数字信号处理实验》指导书

实验一 信号、系统及系统响应

一、实验目的

1、 熟悉连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解。 2、 熟悉时域离散系统的时域特性。 3、 利用卷积方法观察分析系统的时域特性。

4、 掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶变换对连续信号、离散信号及系统

响应进行频域分析。

二、实验内容及步骤

1、 认真复习采样理论、离散信号与系统、线性卷积、序列的傅里叶变换及性质等有关内容，阅读

本实验原理与方法。

2、 了解要使用到的MATLAB命令：

1）基于DTFT离散时间信号分析函数：freqz，real，imag，abs，angle，unwrap。 函数freqz可以用来计算一个以ej?的有理分式形式给出的序列的DTFT值。freqz的形式多样，

常见的有H=freqz(num,den,w),其中num表示序列有理分式DTFT的分子多项式系数，den表示分母多项式系数（均按z的降幂排列），矢量w表示在0～?中给定的一系列频率点集合。freqz函数的其他形式参见帮助。在求出DTFT值后，可以使用函数real, imag, abs和angle分别求出并绘出其实部，虚部，幅度和相位谱。如果需要，还可以用unwrap函数消除相位中的跳变。 2）函数fft(x)可以计算R点序列的R点DFT值；而fft(x,N)则计算R点序列的N点DFT，若R>N，则直接截取R点DFT的前N点，若R

1）产生实验中要用到的下列信号序列： a．采样信号序列：对下面连续信号 xa(t)?Ae?atsin?(0t)u(t)

进行采样，可得到采样序列

xa(n)?xa(nT)?Ae?anTsin(?0nT)u(n)，0?n?50

其中A为幅度因子，a为衰减因子，?0是模拟角频率，T为采样间隔。这些参数都要在实验过程中由键盘输入，产生不同的xa?t?和xa?n?。 b．单位脉冲序列： xb(n)??(n)

1

c．矩形序列： xc(n)?RN(n)，N?10 2）产生系统单位脉冲响应序列。 a． ha(n)?R10(n)

b． hb(n)??(n)?2.5?(n?1)?2.5?(n?2)??(n?3)

3）编写有限长序列线性卷积子程序，用于完成两个给定长度的序列的卷积。也可以直接调用

MATLAB语言中的卷积函数conv。conv用于两个有限长度序列的卷积，它假定两个序列都从n=0开始。调用格式如下： y=conv（x，h）

其中参数x和h是两个已赋值的行向量序列。 4、完成下述实验内容：

1）分析采样序列的特性。产生采样序列xa(n)，A?444.128，a?222，?0?222。 a、 取采样频率fs?1kHz，即T?1ms。观察所采样xa(n)的幅频特性Xe??和xj?a(t)的幅

频特性X?j??在折叠频率处有无明显差别。应当注意，实验中所得频谱是用序列的傅立叶变换公式求得的，所以在频率量度上存在关系：???T。 b、 改变采样频率，fs?300Hz，观察Xe??的变化并做记录。

j?c、 进一步降低采样频率，fs?200Hz，观察频谱混叠是否明显存在，说明原因，并记录

X?ej??。

2） 离散信号、系统和系统响应分析。

a、观察信号xb(n)和系统hb(n)的时域和频域持性；利用线性卷积求信号xb(n)通过系统hb(n)的响应y(n)，比较所求响应y(n)和hb(n)的时域及频域特性，注意它们之间有无差异，绘图说明，并用所学结论解释所得结果。

b、 观察系统ha(n)对信号xa(n)的响应特性。利用线性卷积求系统响应y(n)，并判断y(n)图形

及其非零值序列长度是否与理论结果一致，对xa(n)?ha?n??R10?n?，说出一种定性判断y(n)图形是否正确的方法。调用序列傅立叶变换子程序，求得Ye?j?k?，观察Y?e?j?R的特性曲线，定性判断结果的正确性。改变xa(n)的长度，取N=5，重复该试验。注意参数变化的影响，说明变化前后的差异，并解释所得的结果。

3） 卷积定理验证。将实验2）中的信号换成xa(n)，使a=0.4，?0?2.0734，A=1，T=1，重复

2

实验2）a，打印Ye并绘出Ye卷积定理。

?j?k?曲线；对主程序作简单修改，计算Y?e??X?e?H?e?，

j?kaj?kbj?k?j?k?曲线，与前面直接对进行傅立叶变换所得幅频特性曲线进行比较，验证时域

三、思考题

1、 在分析理想采样序列特性的实验中；采样频率不同时，相应理想采样序列的傅里叶变换谱的数

学频率度量是否都相同?它们所对应的模拟频率是否相同?为什么?

2、 在卷积定理验证的实验中，如果选用不同的频域采样点数M值，例如，选M=10和M=20，分

别做序列的傅里叶变换，求得 Y(ej?k)?Xa(ej?k)Hb(ej?k)，k=0，1，……M-1

所得结果之间有无差异？为什么？

四、实验报告要求

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

2、 按实验步骤附上实验过程中的信号序列、系统单位脉冲响应及系统响应序列的时域和幅频特性

曲线，并对所得结果进行分析和解释。 3、 总结实验中的主要结论。 4、 简要回答思考题。

3

实验二 用FFT作谱分析

一、实验目的

1、 进一步加深DFT算法原理和基本性质的理解（因为FFT只是DFT的一种快速算法，所以FFT

的运算结果必然满足DFT的基本性质）。 2、 熟悉FFT算法原理和FFT子程序的应用。

3、 学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，

以便在实际中正确应用FFT。

二、实验内容及步骤

1、 复习DFT的定义、性质和用DFT作谱分析的有关内容。

2、 复习FFT算法原理与编程思想，并对照DIT—FFT运算流图和程序框图，读懂本实验提供的FFT

子程序。

3、 编制信号产生子程序，产生以下典型信号供谱分析用：

x1(n)?R4(n)

0?n?3?n?1? x2(n)??8?n 4?n?7

?0其它n?0?n?3?4?n?x3(n)??n?3 4?n?7

?0其它n? x4(n)?cos?4n n

x5(n)?sin?8x6(t)?cos8?t?cos16?t?cos20?t

应当注意，如果给出的是连续信号xa(t)，则首先要根据其最高频率确定采样频率fs，以及由频率分辨力确定采样点数N，然后对其进行软件采样（即计算x(n)?xa(nT)，0?n?N?1），产生对应序列x(n)。对信号x6(t)，频率分辨力的选择要以能分辨开三个频率对应的谱线为准则。对周期序列，最好截取周期的整数倍进行谱分析，否则可能产生较大的分析误差，请同学们根据DFT的隐含周期性思考这个问题。 4、完成下述实验内容：

4

1） 对2中所给出的信号逐个进行谱分析。下面给出各个信号的FFT变换区间N以及连续信号x6(t)的采样频率，供实验时参考。 x1(n)：N=8，64 x2(n)，x3(n)：N=8，16

x4(n)，x5(n)：N=8，16

2）令x7(n)?x4(n)?x5(n)，用FFT计算8点和16点离散傅里叶变换。 3）令x8(n)?x4(n)?jx5(n)，用FFT计算8点和16点离散傅里叶变换。

三、思考题

1、在N=8时，x2(n)和x3(n)的幅频特性会相同吗？为什么？N=16呢？ 2、如果周期信号的周期预先不知道，如何使用FFT进行谱分析？

四、实验报告要求

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

2、结合实验中所得给定典型序列幅频特性曲线，与理论结果比较，并分析说明误差产生的原因以及用FFT作谱分析时有关参数的选择方法。 3、总结实验所得主要结论。 4、简要回答思考题。

5

实验三 用双线性变换法设计IIR数字滤波器

一、实验目的

1、熟悉用双线性变换法设计IIR数字滤波器的原理与方法。 2、掌握数字滤波器的计算机仿真方法。

3、通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。

二、实验内容

1、 用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。设计指标参数为：在通带内频率低于

0.2?时，最大衰减小于1dB；在阻带内[0.3?，?]频率区间上，最小衰减大于15dB。 2、 以T=1s为采样间隔，打印出数字滤波器在频率区间[0，?2]上的幅频响应特性曲线。 3、 用所设计的滤波器对实际心电图信号采样序列(在本实验后面给出)进行仿真滤波处理，并分别打

印出滤波前后的心电图信号波形图，观察总结滤波作用与效果。

三、实验步骤

1、 复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR数字滤波器的内容，用双线性变换法设计数字滤波器系统函数H(z)。 教材中满足本实验要求的数字滤波器函数

H?z??3?1?1.2686z?1?0.7051z?2??1?1.010z?1?0.3583z?2??1?1.9044z?1?0.2155z?2?0.007378?1?z?1?5

??Hk?z?k?1A1?2z?1?z?2式中：Hk?z??，k=1，2，3 ?1?21?Bkz?CkzA=0.09036，

B1=1.2686， C1=-0.7051 B2=1.0106， C2=-0.3583 B3=0.9044， C3=-0.2155

根据设计指标，调用MATLAB信号处理工具函数buttord和butter，也可得到H?z?。 由滤波器的函数可见，滤波器H?z?由三个二阶滤波器H1?z?、H2?z?、H3?z?级联组成。 2、 编写滤波器仿真程序，计算H?z?对心电图信号采样序列xa(n) 的响应序列y(n)。

?? 6

yk?n?为第k阶滤波器Hk?z?的输出序列，yk?1?n?为输入序列，根据滤波器H?z?的组成可得

差分方程

yk?n??Ayk?1?n??2Ayk?1?n?1??Ayk?1?n?2??Bkyk?n?1??Ckyk?n?2? 当k?1时，yk?1?n??x?n?。所以H?z?对x?n?的总响序列y?n?可以用顺序迭代算法得到。即依次对k=1，

2，3求解差分方程，最后得到y3?n??y?n?。仿真程序就是实现上述求解差分方程和顺序迭代算法的通用程序，也可直接调用MATLAB filter函数实现仿真。 3、 在通用计算机上运行仿真滤波程序，完成实验内容（2）和（3）。

四、实验报告要求

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

j?2、 由所打印的H(e)特性曲线及设计过程简述双线性变换法的特点。

五、心电图信号采样序列x?n?：

人体心电图信号在测量过程中往往受到工业高频干扰，所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。下面给出一实际心电图信号采样序列样本x(n)，其中存在高频干扰。在实验中，以x(n)作为输入序列，滤除其中的干扰成分。

x(n) ={-4, -2, 0, -4, -6, -4, -2, -4, -6, -6,

-4, -4, -6, -6, -2, 6, 12, 8, 0, -16, -38, -60, -84, -90, -66, -32, -4, -2, -4, 8, 12, 12, 10, 6, 6, 6, 4, 0, 0, 0, 0, 0, -2, -4, 0, 0, 0, -2, -2, 0, 0, -2, -2, -2, -2, 0}

7

实验四 用窗函数法设计FIR数字滤波器

一、实验目的：

1、掌握用窗函数法设计FIR数字滤波器的原理和方法。 2、熟悉线性相位FIR数字滤波器特性。 3、了解各种窗函数对滤波特性的影响。

二、实验内容及步骤：

1、复习用窗函数法设计FIR数字滤波器一节内容。

一般，设计线性相位FIR数字滤波器采用窗函数法或频率抽样法，本实验采用窗函数法，采用矩形窗。下面简要介绍窗函数法设计FIR理论知识，更详细的相关知识请参考教材。

如果所希望的滤波器的理想频率响应函数为Hd(ejw)，如理想的低通，由信号系统的知识知道，在时域系统的冲激响应hd(n)将是无限长的，如图1、图2所示。 Hd(w) -wc wc

图2

图1

若时域响应是无限长的，则不可能实现，因此需要对其截断，即设计一个FIR滤波器频率响应

H(e)??h(n)e?jwn来逼近Hd(ejw)，即用一个窗函数w(n)来截断hd(n)，如式3所示：

jwn?0N?1 h(n)?hd(n)w(n)（式3）。

最简单的截断方法是矩形窗，实际操作中，直接取hd(n)的主要数据即可。

h(n)作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数为：

H(e)??h(n)e?jwn（式4）

jwn?0N?1令z?e，则

，式中，N为所选窗函数w(n)的长度。 H(z)??h(n)z?n（式5）

n?0N?1jw 8

如果要求线性相位特性，h(n)还必须满足：

h(n)??h(N?1?n)

（式6），根据式6中的正、负和长度N的奇偶性又将线性相位FIR

滤波器分成四类。要根据所设计的滤波器特性正确选择其中一类。例如：要设计线性相位低通特性，可选择h(n)?h(N?1?n)类。

有关窗函数的种类及特性，请参照教材。

2、编写程序

1）编写能产生矩形窗、升余弦窗、改进升余弦窗和二阶升余弦窗的窗函数子程序，也可调用MATLAB相关函数。

2）编写低通FIR数字滤波器主程序。 3）上机实验内容。

A、用升余弦窗设计一线性相位低通FIR数字滤波器，截止频率?c??4rad。窗口长度N＝15,33。

要求在两种窗口长度情况下，分别求出h(n)，打印出相应的幅频特性和相频特性曲线，观察3dB带宽和20dB带宽。总结窗口长度N对滤波特性的影响。 B、N＝33,?c＝

?，用四种窗函数设计线性相位低通滤波器。绘制相应的幅频特性曲线，观4察3dB和20dB带宽以及阻带最小衰减，比较四种窗函数对滤波器特性的影响。

三、思考题：

1、 如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通

滤波器?写出设计步骤。

2、 如果要求用窗函数法设计带通滤波器，且给定上、下边带截止频率为?1和?2，试求理想带通

的单位脉冲响应hd(n)。

四、实验报告要求：

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

2、 按照实验步骤及要求，比较各种情况下的滤波性能，说明窗口长度N和窗函数类型对滤波特性

的影响。

3、 总结用窗函数法设计FIR滤波器的主要特点。 4、 简要回答思考题。

9

参考：用于数字信号处理的MATLAB操作或者函数

1、向量的生成：利用冒号（：）生成向量。 1）X=j:k (j

生成X=[j,j+1,j+2,…,k-1,k] 2）X=j:I:k

如果I>0 且jk,则生成向量

X=[j,j+I,j+2I,…,k-I,k] 如果I<0 且j>k,则生成向量 X=[j,j+I,j+2I,…,k]

例： X1=1：5 X1=[1 2 3 4 5]

X2=1：0.5：3 X2=[1.000 1.500 2.000 2.500 3.000] X3=5：-1：1 X3=[5 4 3 2 1]

2、矩阵的生成：

zeros生成全0阵

B= zeros(n) 生成n*n的全0矩阵 B= aeros(m,n) 生成m*n的全0矩阵

B= zeros(size(A)) 生成与与矩阵A大小相同的全0矩阵 注：m,n必须是一个非负数 ones生成全1阵

rand生成均匀分布的随机阵 3、矩阵的算术运算 1）加法和减法

对于同维矩阵指令为：A+B、A-B；

对于矩阵和标量（一个数）的加减运算，指令为：A+3、A-9 2） 乘法和除法运算

A*B 是数学中的矩阵乘法，遵循矩阵乘法规则； A．*B 是同维矩阵对应位置元素做乘法； B=inv（A）是求矩阵的逆；

A/B 是数学中的矩阵除法，遵循矩阵除法规则； A．/B 是同维矩阵对应位置元素相除；

A' 表示矩阵的转置运算；

10

4、数组函数

下面列举一些基本函数，他们的用法和格式都相同。 sin(A),cos(A),exp(A),log(A)(相当于ln)

sqrt(A)开平方 abs(A)求模 real(A)求实部 imag(A)求虚部

式中A可以是标量也可以是矩阵 例： 利用等差向量产生一个正弦值向量 t=0:0.1:10 A=sin(t) plot(A)

这时候即可看到一个绘有正弦曲线的窗口弹出。

另：每条语句后面加“；”表示不要显示当前语句的执行结果； 不加“；”表示要显示当前语句的执行结果。

5、绘图（二维绘图）

利用前例

如果要使向量的横纵坐标一一对应，则应写为plot(t,A)。 如不要横坐标对应时，可只写plot(A)。

如果A是一个矩阵不是向量，则在同一窗口中可绘出与矩阵行数相同的曲线且颜色不同。 同时也可指定曲线的颜色和格式。 如plot(t,A,’b.’)

则原来的连续曲线就变成了蓝色的点线。 下面给出常用选项。

r红 g绿 b蓝 y 黄 m 洋红 c青 w白 k黑 ━实线 ┅┅虚线

*用星号标数据点；.用点标；o（字母“o”）用圆圈标；x（字母“x”）用叉标；-.点划线 下面给出窗口分割语句:subplot(x,y,n)表示将一个窗口分割成x行，y列 ，共x*y个窗口， 取第n个窗口为当前窗口。

例subplot(3,2,2) 分割3*2个窗口，取第二个窗口，如果下面有绘图语句，就表示要在第二个窗口中绘图。 6、联机帮助

1)help —准确指令，如 help elfun help exp help ftt

2)lookfor —有关指令，如 lookfor fourier

11

本实验可能用到的一些序列和函数

1）单位冲击序列??n?N?：

function [x,n]=impseq(n0,n1,n2) n=[n1:n2];

x=[(n-n0)==0];利用逻辑运算实现在n=N处为1。如：impseq(5,0,10) 2）单位阶跃序列u?n?N）?：

function [x,n]=stepseq(n0,n1,n2) n=[n1:n2]

x=[(n-n0)>=0]; 利用逻辑运算实现在n>=N处为1。如：stepseq(5,0,10) 3）矩形序列：

function [x,n]=stepseq(n0,n1,n2) n=[n1:n2]

x=[(n-n0)>=0] stepseq(0,0,9)实现了R10?n?

也可以用

function [x,n]=squseq(M，N,n1,n2) n=[n1:n2]

x=[(n-M)>=0&(n-N)<0] 利用逻辑运算实现在N>n>=M处为1 4) 正弦序列

x?n??sin?n*??

x(n)=sin(n*pi) x(t)=sin(2*t)

5) 也可调用库函数产生一些信号如sinc 6）直接输入序列或者函数表达式 x(t)=sin(2*t)

x(n)=[1,1,1,1,1,1]

掌握二维绘图方法 1）plot —连续曲线 t=0:pi/50pi; y=sin(t); plot(t,y,)

hold on保留前面坐标系中已经存在的图形 plot(t,y,s1,s2,s3) s1,s2,s3表示线型、颜色等 2）stem离散图形 [n,y]=squseq(4,4,0,10); stem(n,y,’o’)

3)一个窗口中绘制多个图形 subplot(2,3,1) plot(x)

subplot(2,3,2) plot(y)

12

函数调用

常用的函数：fft，abs，angle及滤波器设计中的函数见下面滤波器的设计中 数字巴特沃兹滤波器的设计方法 方法一

%获得模拟滤波器的性能指标，边界频率进行预畸变校正，按公式??Wp=0.2*pi Ws=0.3*pi;

Rp=1;%单位DB Rs=15; %单位DB Fs=0.02*pi; %单位Hz %模拟滤波器设计

[n,Wn]=buttord(Wp,Ws,Rp,Rs);%获得阶数和通带截止频率 [z,p,k]=buttap(n);%设计模拟低通滤波器原型

[b0,a0]=zp2tf(z,p,k);%转化为传递函数表示的滤波器

[b,a]=lp2lp(b0,a0,Wn);%将滤波器原形转化为要求设计的低通滤波器 %使用双线性变换法得到设计的数字滤波器 [bz,az]=bilinear(b,a,Fs);

[h,w]=freqz(bz,az);%计算频响 绘频率响应曲线 plot(h); 方法二

1、 数字滤波器的性能指标归一化 Wp=0.2 Ws=0.3

Rp=1;%单位DB Rs=15; %单位DB

2、 获得滤波器的阶数和通带截止频率 [n,Wn]=buttord(Wp,Ws,Rp,Rs) 设计巴特沃兹低通数字滤波器 [b,a]＝butter(n,Wn) freqz=(b,a)

方法三见本实验附录程序

A=[0.09036,2*0.09036,0.09036]; B1=[1,-1.2686,0.7051]; B2=[1,-1.0106,0.3583]; B3=[1,-0.9044,0.2155]; [H1,w]=frez(A,B1,100); [H2,w]=frez(A,B2,100); [H3,w]=frez(A,B3,100); H=H1.*H2.*H3; mag=abs(H);

db=20log((mag+eps)/max(mag)); subplot(2,2,1) plot(w/pi,db)

13

2???tan??进行 T?2?DFT的源程序 functionX=dft(x) N=length(x); n=0:N-1; k=0:N-1;

WN=exp(-j*2*pi/N);%求得WN因子 Nk=n`*k;

WNnk=WN.^nk;%求得DFT运算矩阵，注意向量乘为点乘 X=x*WNnk;%利用矩阵运算计算DFT

14

本文来源：https://www.bwwdw.com/article/qebo.html