数字信号处理实验（上实验课时请同时带教材和实验指导书）
南昌航空大学信息工程学院
2014年3月
目录
实验一熟悉MatLab (1)
一、实验目的 (1)
二、Matlab基本练习 (1)
（一）MA TLAB功能简介 (1)
（二）MatLab中的向量及矩阵表示练习 (1)
（三）MatLab基本绘图 (3)
（四）MatLab帮助的使用 (5)
（五）m文件的建立和调试 (6)
三、信号处理相关综合练习 (7)
四、实验内容 (7)
五、实验报告内容及要求 (7)
实验一参考例程 (7)
实验二时域采样定理 (8)
一、实验目的 (8)
二、实验原理 (8)
三、实验用MATLAB函数介绍 (9)
四、实验内容和步骤 (9)
五、实验报告要求 (10)
参考程序 (10)
实验三用DFT(FFT)对信号进行谱分析 (12)
一、实验目的 (12)
二、实验原理 (12)
三、实验内容及步骤 (14)
四、实验报告要求 (15)
参考程序 (15)
实验四用双线性变换法设计IIR数字滤波器 (21)
一、实验目的 (21)
二、实验原理 (21)
三、实验内容 (21)
四、实验步骤 (22)
实验四参考例程 (22)
实验五用窗函数法设计FIR数字滤波器 (23)
一、实验目的 (23)
二、FIR低通数字滤波器的设计原理 (23)
三、实验内容 (24)
四、实验步骤 (24)
实验五参考例程 (25)
实验一熟悉MatLab
一、实验目的
1、熟悉MATLAB软件的集成开发环境，学会利用MATLAB编程及获得帮助的方法。
2、掌握MATLAB中有关矩阵和向量的表示和运算方法。
3、学会利用MATLAB的绘图功能。
4、学会对MA TLAB程序进行调试
二、Matlab基本练习
（一）MATLAB功能简介
在科学研究和工程应用中，往往要进行大量的数学计算，这些运算一般来说难以用
手工精确和快捷地进行，而要借助计算机编制相应的程序做近似计算。
Matlab目前已经成为国际上流行的仿真软件之一，除了传统的交互式编程之外，还提供了丰富可靠的矩阵运算、图形绘制、数据处理、方便的Windows编程等便利工具。
Matlab的应用领域涉及到各个方面，而且形成了专门的工具箱，如控制系统工具箱（control systems toolbox）、系统辨识工具箱(system identification toolbox)、鲁棒控制工具箱(robust control toolbox)、神经网络工具箱（neural network toolbox）、最优化工具箱（optimization toolbox）、信号处理工具箱(signal processing toolbox)、小波分析工具箱(wavelet toolbox)以及仿真环境SIMULINK等。
由于Matlab编程运算方式不像学习其它高级语言（如Basic、Fortran和C等）那样难于掌握。Matlab语言具有编程效率高、易学易懂、使用方便、绘图方便等一系列优点。
（二）MatLab中的向量及矩阵表示练习
命令窗口是MA TLAB的主要交互窗口，用于输入命令并显示除图形以外的所有执行结果。
MATLAB命令窗口中的“>>”为命令提示符，表示MA TLAB正在处于准备状态。在命令提示符后键入命令并按下回车键后，MA TLAB就会解释执行所输入的命令，并在命令后面给出计算结果。
（1）矩阵输入：
在MatLab环境下，矩阵输入的方式很直观，如矩阵
A的输入方式为：
A=[1,2,3;4,5,6;7,8,9]，其中同一行的内容用逗号或空格隔开，而采用分号表示换行。
按这种方式输入矩阵后，按回车键，矩阵A的内容将在屏幕上显示出来。
A=[1,2,3;4,5,6;7,8,9]
A =
1 2 3
4 5 6
7 8 9
如果在上式所示语句的后面加一个分号，则其结果便不会在屏幕上显示。如
A=[1,2,3;4,5,6;7,8,9]；
（2）向量输入：
1）向量的直接输入，例如输入下面两条命令
V1=[1 2 3 4]
V1=[1;2;3;4]
可以简单的输出下面的行向量和列向量
2）等差向量的输入：向量除了可以直接输入外，还可以利用“：”运算符获得等差向量。
如输入1：0.1：2则会产生从1开始一直到2公差为0.1的向量，即[1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000 1.8000 1.9000 2.0000]。
练习：
x=1:0.1:2
x =
Columns 1 through 8
1.0000 1.1000 1.2000 1.3000 1.4000 1.5000 1.6000 1.7000
Columns 9 through 11
1.8000 1.9000
2.0000
当公差为1时可以省略，如1：5；产生的向量为[1 2 3 4 5]；
n=1:5
n =
1 2 3 4 5
（3）零向量、零矩阵、1向量、1矩阵
与向量生成有关的简单函数还有zeros(m,n)产生m×n维0矩阵;ones(m,n)产生m×n维1矩阵。
x=ones(2,3)
x =
1 1 1
1 1 1
（4）数组
在MATLAB中可以采用如下的形式引用数组元素。A(:,2)表示A矩阵的第2列元素；A(2,:)表示A矩阵中的第2行元素；A(1,2)表示矩阵中的第1行第2列的元素。（注：在MA TLAB 中矩阵元素的下标从1开始）。
如输入：A=[1 2 3;4 5 6;7 8 9];
A(2,:)
输出：ans=
4 5 6
输入：A(1,2)
输出：ans=
2
（三）MatLab基本绘图
绘图涉及到两个方面，指定自变量的范围，选择绘图命令及绘图格式。
（1）stem函数
“stem”绘制针状图。“stem”适合绘制离散序列的图形。例如，下面的代码可以绘制正弦函数的图形：
n=1:50;%定义序列的长度是50
t=2*pi*n/50;%定义时间变量，自变量范围为[0,2*pi]
xn=sin(t);%计算正弦函数值
stem(xn);%画针状图
得到的函数图形如图1.1所示：
图1.1针状图示例
（2）plot函数
“plot”函数绘制图形，其调用格式一般为plot(x,y)，其中x和y分别是欲画图像在X轴和Y轴的数据向量。点及其连接类型有多种选项（查帮助可得）。点数足够密时，成为光滑曲线。如输入以下内容，可以绘制正弦函数等的图形，并对图形处理、加注释等：close all; %关闭所有已经打开的图形视窗
x = -pi:0.1:pi; %指定自变量的范围
y = sin(x);
plot(x,y) ; %绘制正弦曲线，如图1.2
figure; %打开新的图形窗口。去掉该语句，新的图形将覆盖旧的图形
plot(x, sin(x), x, cos(x)); %同时绘制两条曲线,图1.3
axis([-3.2,3.2, -1.2, 1.2]); %axis([xmin,xmax,ymin,ymax])函数来调整图轴的范围xlabel('Input Value'); %x轴注解
ylabel('Function Value'); %y轴注解
title('Two Trigonometric Functions'); %图形标题
legend('y = sin(x)','y = cos(x)'); %图形注解
grid on; % 显示格线，最后图形见图1.4
图1.2 正弦曲线
图1.3 绘制两条正弦曲线
图1.4 加注释等之后的图1.3
(3)函数“subplot”，该函数用来同时画出数个小图形于同一个视窗之中，例如：
close all; %关闭所有已经打开的图形视窗
x = -pi:0.1:pi; %指定自变量的范围
subplot(2,2,1); plot(x, sin(x)); %一个大图中分四个小图，两行、两列，此为1号位置图subplot(2,2,2); plot(x, cos(x)); %一个大图中分四个小图，两行、两列，此为2号位置图subplot(2,2,3); plot(x, sinh(x));
subplot(2,2,4); plot(x, cosh(x));
图1.4 subplot示例
（四）MatLab帮助的使用
1．MA TLAB帮助文档：主窗口下按F1打开帮助文档，查找、搜索相关内容
Contents一栏中的MATLAB下的Getting Started中有基本的MA TLAB命令介绍
signal processing toolbox中有本课程所用的各种信号处理相关命令
2．快捷键获得帮助
图1.5 获取命令帮助
3．help 命令名获取帮助，如：
help plot
（五）m文件的建立和调试
（1）建立新的M文件
1) 菜单操作。从MATLAB主窗口的File菜单中选择New菜单项，再选择M-file命令，屏幕上将现MA TLAB 文本编辑器窗口,输入M文件的内容并存盘。
2) 命令按钮操作。单击MA TLAB主窗口工具栏上的New M-File命令按钮，启动MATLAB文本编辑器.
（2）调试：在程序当中设置断点，或单步执行，可观察每步的运行结果。
图1.6 m文件的调试
?F5：运行到断点
?F10：单步运行
?设置、取消断点：鼠标左键在左边行号处点击
?工作空间是MA TLAB用于存储各种变量和结果的内存空间。调试时，可在工作空间
对变量进行观察、编辑、保存和删除。简单变量的值可以用鼠标放在m文件的相应变量上进行显示。
三、信号处理相关综合练习
（1）建立m文件，执行卷积运算，进行设置断点、调试练习，并绘针状图
用MATLAB计算序列{-2 0 1 –1 3}和序列{1 2 0 -1}的离散卷积。在实验报告中给出卷积结果和图形结果
函数conv功能简介：计算离散卷积。
调用格式：y n=conv(hn, xn)计算向量hn和xn的卷积，卷积结果向量为yn。如果用N 和M分别表示hn和xn的长度，则yn的长度为N+M-1。
（2）实现单位抽样序列δ(n)和矩形序列R10(n-2)。序列起始序号为-5, 结束序号为20。要求绘出杆状图。
四、实验内容
1、熟悉MATLAB软件的集成开发环境，学会利用MATLAB编程及获得帮助的方法。
2、做内容二中的练习。
3、做内容三的练习，并在报告中给出实验程序和实验结果
五、实验报告内容及要求
1、简述实验目的。
2、内容三的题目、程序及执行结果，要求形成word文档并打印（matlab图形中用edit －copy figure菜单可将图形拷贝到word中）。
实验一参考例程
三、信号处理相关综合练习
（1）卷积运算
function experiment11 % 函数名，一般与文件名相同
close all; %关闭所有已经打开的图形视窗
a=[-2 0 1 -1 3];
b=[1 2 0 -1];
c=conv(a,b) % 卷积函数，求得的结果存放在数组c中
M=length(c)-1; % 计算长度，得出绘图的终止点
n=0:1:M;
stem(n,c);
xlabel('n'); ylabel('幅度'); （2）单位抽样序列：
close all n1=-5 n2=20 n0=0 n=[n1:n2]
x = zeros(1,n2-n1+1) x(n0-n1+1)=1 stem(n,x) R 10(n-2)的实现：
close all n0=-5 n1=2 n2=11 n3=20 n=[n0:n3]
x = zeros(1,n3-n0+1) x(n1-n0+1:n2-n0+1)=1 stem(n,x)
实验二 时域采样定理
一、 实验目的
熟悉并加深对采样定理的理解，了解采样信号的频谱和模拟信号频谱之间的关系。
二、 实验原理
模拟信号经过理想采样，形成采样信号。采样信号的频谱和模拟信号频谱之间的关系如下：
对模拟信号x a (t )进行理想采样的公式如下式：
对上式进行傅立叶变换，得到：
1?()(())a a s m X j X j m T ∞
=-∞Ω=Ω-Ω∑
?()()()a a m x
t x t t nT δ∞
=-∞
=-∑
上式的右边就是序列的傅立叶变换X (e j ω)，即
三、 实验用MATLAB 函数介绍
1. fft
功能： 一维快速傅立叶变换(FFT)。
Xk=fft(xn, N)： 采用FFT 算法计算序列向量xn 的N 点DFT 。缺省N 时, fft 函数自动按xn 的长度计算xn 的DFT 。当N 为2的整数次幂时，fft 按基2算法计算，否则用混合基算法。
2. ifft
功能： 一维快速逆傅立叶变换(IFFT)
。
调用格式: 与fft 相同。 四、
实验内容和步骤
1. 给定模拟信号如下：
x a (t )=Ae -αt sin(Ω0t )u (t )
x a (t )的幅频特性曲线
2. 按照选定的采样频率对模拟信号进行采样，得到时域离散信号x (n ):
x (n )=x a (nT )=Ae -αnT sin(Ω0nT )u (nT )
这里给定采样频率如下： f s =1 kHz ，300 Hz ，200 Hz 。分别用这些采样频率形成时域离散信号, 按顺序分别用x 1(n )、 x 2(n )、 x 3(n )表示。选择观测时间T p =50 ms 。
3. 计算x (n )的傅立叶变换X (e j ω
):
??()()[()()]j t a a j t a m X j x
t e dt x t t nT e dt
δ∞
-Ω-∞∞
∞
-Ω-∞
=-∞
Ω==-?∑??()()|j a t
X j X e ωω=ΩΩ=10.8
0.60.40.200
100
200300
400
500
f / Hz
|X a (j f )|
1
00
()[()]sin()i n j at j n
i n X e FT x n Ae nT e ω
ω---===Ω∑123111(,,)T s T s T s =
==
式中, i =1, 2, 3，分别对应三种采样频率的情况 采样点数用下式计算：
因ω是连续变量。为用计算机进行数值计算，改用下式计算：
式中，ωk = (2πMk, k=0, 1, 2, 3, …, M-1； M=64。可以调用MATLAB 函数fft
计算
4. 打印三种采样频率的幅度曲线|X (e j ωk )|~ωk , k=0, 1, 2, 3, …, M-1，M=64。
五、 实验报告要求
1. 简述实验原理。
2.
针对三种采样频率情况，进行分析讨论。
3.
参考程序
function experiment2
%采样定理演示：频谱混叠问题 close all;%
% 总的参数 A=444.128; a=50*pi*2^0.5; W0=50*pi*2^0.5
Tp=0.05 %此处波形趋于零，信号时间截取到此为止，50毫秒
% 模拟信号频谱的离散近似计算，即信号的完整频谱
Ts = 0.00005;%采样周期，取较小间隔，获取较大采样频率，近似实现连续信号频谱的效果,fs=1/Ts=20000
fs=1/Ts; t = 0:Ts:Tp;
xn = A*exp(-a*t).*sin(W0*t); N=Tp/Ts+1 k=0:N-1
f=fs/N*k %f=fs/N*k ，2*pi/N*k=2*pi*f/fs,对应各个k 值的实际频率值,换算为实际频率 Xk=fft(xn);
subplot(4,2,1);plot(t,xn);
p i i
T n T =
100
()[()]sin()i k
i k n j anT j n
M i n X e
DFT x n Ae nT e ωω---===Ω∑
xlabel('模拟信号时间单位：秒'); ylabel('信号')
subplot(4,2,2);plot(f,abs(Xk));
xlabel('频率：Hz'); ylabel('模拟信号频谱')
xlim([0,500]);%500Hz频率频谱接近零，因此显示到此处
ylim([0,max(abs(Xk))+0.05*max(abs(Xk))]);%y轴显示比最高值稍大一点，便于看清顶点
% 1000Hz采样
Ts = 1/1000;%fs=1000Hz,采样点数较少
fs=1/Ts;
t = 0:Ts:Tp;
xn = A*exp(-a*t).*sin(W0*t);
M=Tp/Ts+1;N=64;%计算采样点数，小于64
k=0:N-1 %计算64点fft
f=fs/N*k %f=fs/N*k，2*pi/N*k=2*pi*f/fs,对应各个k值的实际频率值,换算为实际频率Xk=fft(xn,64);
subplot(4,2,3);stem(t,xn);
xlabel('采样信号时间单位：秒'); ylabel('信号')
subplot(4,2,4);plot(f,abs(Xk));
xlabel('频率：Hz'); ylabel('1000Hz采样频谱')
xlim([0,500]);%显示到此处，显示了半个周期
ylim([0,max(abs(Xk))+0.05*max(abs(Xk))]);%y轴显示比最高值稍大一点，便于看清顶点
% 300Hz采样
Ts = 1/300;%fs=1000Hz,采样点数较少
fs=1/Ts;
t = 0:Ts:Tp;
xn = A*exp(-a*t).*sin(W0*t);
M=Tp/Ts+1;N=64 ;%计算采样点数，小于64
k=0:N-1 %计算64点fft
f=fs/N*k %f=fs/N*k，2*pi/k=2*pi*f/fs,对应各个k值的实际频率值,换算为实际频率
Xk=fft(xn,64);
subplot(4,2,5);stem(t,xn);
xlabel('采样信号时间单位：秒'); ylabel('信号')
subplot(4,2,6);plot(f,abs(Xk));
xlabel('频率：Hz'); ylabel('300Hz采样频谱')
xlim([0,500]);%300Hz之后未计算
ylim([0,max(abs(Xk))+0.05*max(abs(Xk))]);%y轴显示比最高值稍大一点，便于看清顶点
% 200Hz采样
Ts = 1/200;%fs=1000Hz,采样点数较少
fs=1/Ts;
t = 0:Ts:Tp;
xn = A*exp(-a*t).*sin(W0*t);
M=Tp/Ts+1;N=64 ;%计算采样点数，小于64
k=0:N-1 %计算64点fft
f=fs/N*k %f=fs/N*k ，对应各个k 值的实际频率值,换算为实际频率 Xk=fft(xn,64);
subplot(4,2,7);stem(t,xn);
xlabel('采样信号 时间单位：秒'); ylabel('信号')
subplot(4,2,8);plot(f,abs(Xk)); %计算频率响应的模的值并绘图。 xlabel('频率：Hz'); ylabel('200Hz 采样频谱') xlim([0,500]);%200Hz 之后未计算
ylim([0,max(abs(Xk))+0.05*max(abs(Xk))]);%y 轴显示比最高值稍大一点，便于看清顶点
实验三 用DFT(FFT)对信号进行谱分析
一、实验目的
1、在理论学习的基础上，通过本次实验加深对DFT 的理解。
2、熟悉应用FFT 对典型信号进行频谱分析的方法。
3、了解应用FFT 进行信号频谱分析过程中可能出现的各种误差，以便在实际中正确应用FFT 。
二、实验原理
在运用DFT 进行频谱分析的时候可能会产生三种误差，现分析如下： （一）截断效应
实际中的信号序列往往很长，甚至是无限长序列。为了方便，我们往往只取实际序列的一部分来近似它们。这种截短等价于给原信号序列乘以一个矩形窗函数。根据卷积定理，最终信号的频谱等于原信号的谱和矩形窗的谱的卷积，从而造成谱线加宽或称为频谱泄漏。矩形窗时间取得越长，矩形窗的频谱变窄，由截断引起的效应会减小。
例如50 Hz 正弦波xa (t )=sin(2π·50t)，它的幅度曲线是线状谱，如图3.1(a)所示。如果将它截取0.09s 的一段，相当于将它乘一长度为0.09 s 矩形窗函数，即xa (t )RTp (t)，Tp =0.09s,该信号的谱等于原信号的谱和矩形窗的谱的卷积，如图1（b ）所示。矩形窗长度扩大Tp =0.18s,后，频谱泄漏会变小，如图1（c ）。
1
0.50
－250
－200
－150
－100
－50
050
100
150
200
250
幅度
f / Hz
(a )
1
0.50幅度
－250
－200
－150
－100
－50
050
100
150
200
250
f / Hz
(b )
图 3.1 用DFT 对正弦波进行谱分析
(a)50 Hz 正弦波的幅频曲线；
(b) 50 Hz 正弦波加窗后的幅频曲线(T p=0.09 s); (c) 50 Hz 正弦波加窗后的幅频曲线(T p=0.18 s)
同时，由于频谱泄漏，还会造成靠得很近的两个谱峰混淆为一个谱峰，或是强的谱线的旁瓣掩盖弱的谱线，称为谱间干扰，导致频谱分辨率降低。
矩形窗时间取得越长，矩形窗的频谱变窄，由截断引起的效应会减小。泄漏和谱间干扰都会小。
（二）栅栏效应
因为DFT 是Z 变换在单位圆上的等间隔采样，所以它的频谱不是—个连续函数。这样就产生了栅栏效应。用DFT 来观看频谱就好像通过—个栅栏来观看一幅景象，看到的只是信号在一些离散点上的频谱。这样，如果有一些频谱的值处在离散点之间，就不会被计算，相当于被“栅栏”挡住，不能被我们观察到。减小栅栏效应的—个方法是在原序列的末端补—些零值，从而改变DFT 的点数。这种方法的实质是人为地改变了对真实频谱采样的点数和位置，相当于搬动了“栅栏”的位置，从而使得原来被挡住的—些频谱的峰点或谷点显露出来。但补零对改善谱间干扰无作用。 （三）混叠现象
当采样速率不满足Nyquist 定理，即采样频率小于两倍信号(这里指的是实信号)的最高频率时，由于周期性延拓必然造成频率混叠。所以，在利用DFT 分析连续信号频谱的时侯，必须注意这一问题。实验二的时域采样定理实验反映了混叠现象的特点。
在实际应用中要尽可能减小和消除DFT 用于信号的频谱分析时这些误差的影响。 （四）频谱分析时种因素的综合考虑
在实际用FFT 对模拟信号进行谱分析时，首先要把模拟信号转换成数字信号，转换时要求知道模拟信号的最高截止频率，以便选择满足采样定理的采样频率，避免混叠效应。一般选择采样频率是模拟信号中最高频率的3~4倍。如果模拟信号不是严格的带限信号，会因为频谱混叠现象引起谱分析的误差，这种情况下可以预先将模拟信号进行预滤，或者尽量将采样频率取高一些。
除选择采样频率外，还需要确定分辨率F ，并由此进一步确定观测时间，避免栅栏效应。同时，取样时间加长，可以减小泄漏，从而使得相近频率的谱间干扰减小。
最小的观测时间T pmin 和分辨率成倒数关系，一般用教材(6.7.8)式确定。最小的采样点数用教材(6.7.12)式确定。要求选择的采样点数和观测时间大于它的最小值。
f smin =2 f c
用FFT 作谱分析时，一般取FFT 的点数服从2的整数幂，这一点在上面选择采样点数时可以考虑满足，有时可以通过在序列尾部加0完成。
如果要进行谱分析的模拟信号是周期信号，最好选择观测时间是信号周期的整数倍。如果不知道信号的周期，要尽量选择观测时间长一些，以减少截断效应的影响。举一个极端的
min min 2s C f f
N F F
==min 1p T F
=
例子，一个周期性正弦波，如果所取观察时间太短，例如取小于一个周期，它的波形和正弦波相差太大，肯定误差很大，但如果取得长一些，即使不是周期的倍数，这种截断效应也会小一些。
对一个或数个孤立的正弦信号进行频谱分析，或者说连续频谱中存在一些突出尖锐的峰，则栅栏效应的影响较大，且有一些另外的特性需要考虑。在满足采样频率大于3-4倍的信号频率的条件下，取信号周期的整倍数的时间长度进行谱分析，可以避开泄漏。避开泄漏的原理，其推导请参见胡广书编《数字信号处理》（1997年版）“3.8 关于正弦信号抽样的讨论”。从图形上面分析，可参考图3.1。
回顾相关知识，矩形窗的FFT的主峰及零点在2kπ/N。观察图3.1,单个频率的信号截断后，其频谱是其本来的线谱同矩形窗卷积的结果，频谱以线谱为核心向两边产生泄漏，主峰。做FFT时，频率点的取值是：ωk＝2πf k/f s＝2πk/N，相对应的模拟频率f k =kf s/N。假如信号的频率为f x, 那么如果f x能够和某个f k重合，则FFT的频谱就可以避开截断效应。对正弦信号f x能够和某个f k重合的条件是取信号的整周期进行分析，FFT的k值取样点f k =kf s/N，刚好位于图3.1所示的谱线主峰和周围的零点上（矩形窗的频谱中，零点之间的间隔相等），从而只显示了主峰，泄漏未显示。但当信号的频率未知时，或是存在多条谱线时，这种巧合无法事先推断，比较可靠的方法还是将信号的时间取长一些，长到一定程度之后，再增加时间长度，频谱基本接近。实验中请仔细观察现象。
三、实验内容及步骤
1.初步练习：用信号x1(t)=cos20πt观察频谱泄漏现象，栅栏效应及正弦波抽样时的规律。信号频率为f =10Hz，周期为0.1秒。令采样频率为fs=100Hz，采样周期Ts=0.01s，满足采样定理要求及一般正弦信号抽样习惯。
令截取长度为N=68点，不满足整周期，观察频谱泄漏。根据f k =kf s/N计算可得，模拟信号x1(t)的频率没有与任意一个f k值相对应，而是处在k=6和k=7对应的频率的中间。所以假如没有频谱泄漏，N=68时x1(t)的频谱中不能观察到任何谱线。但由于存在泄漏，观察k=7附近，可以观察到泄漏出的主峰中心线附近的频谱的谱线。
令截取长度为N=100点，共取10个整周期，根据f k =kf s/N计算可得，模拟信号x1(t)对应的k值为10,在k=10处观察到谱线。周围未观察到频谱泄漏。这时候因为截断的关系，泄漏依旧存在，其实质在于泄漏处的不为零的频谱都处在k值的中间，被栅栏挡住未显示（FFT只计算并显示了整数k对应的频谱点）。整数k值刚好处在主谱线以及泄漏频谱的零点上。主谱线被显示，泄漏的非零点未显示。
假如抽样点较少，对于正弦函数，抽样有没有位于整周期处，则由于栅栏效应，k值实际取到的位置距离信号的实际频谱的位置就会太远，不能代表信号的实际频率。因此信号中若含有单个的正弦成分，在不知道正弦成分的准确周期的情况下，抽样点数要增加，直到前后两次做出的频谱接近到一定精度。
由本实验应该看出，单个频率的信号，由于截断效应的影响，向两边泄漏，但存在一个主峰。如果k值抽样对应的频率刚好取在主峰未知，就能避开栅栏效应的影响，较好地反映信号包含的实际频率。
2．进一步练习观察：由包含两个频率的信号x2(t)=cos20πt＋2cos55πt。分别取N=100,512,1024。可以看出，只有取到较大的值，即时间取长之后，才能使两个信号均有效避免栅栏效应的过度影响，反映两个信号的真实相对幅度大小（2倍关系）。如果取样时间过短，那么k值抽样时，就不一定都会抽取在主峰或其附近，造成频谱分析的误差。
3．较复杂的练习：用信号x3(t)=cos0.48πt＋cos0.52πt，其频率分别为0.24Hz, 0.26Hz。
采样周期选1秒，满足采样定理要求。然后分别取（1）N=10；（2）N=10然后补90个零。（3）N=100,101,102;(4) N=512,513;(5)N=20000,20003,20005。分别观察不同N值下的频谱。同时结合参考程序的注释，分析截断、栅栏效应对于频谱的影响。可以看出，不同的N值时，泄漏情况大小不一样，同时，栅栏的位置也不一样，因而造成了频谱的情况不一样。只有N值较大时（本题对应取值时间较长），才能有效避免截断、栅栏效应的影响，但计算代价增加了。
4．给定取样后的信号数据文件experiment34.txt，进行频谱分析。已知采样频率为2Hz。按取样点数从128开始，选取不同的采样点分析6次以上，指出所包含的主要频率及相对幅度。
四、实验报告要求
1.在实验报告中简述实验目的和实验原理要点。
2.在实验报告中附上实验结果，分析所得到的结果。附上所设计程序。
3.总结实验中的主要结论。
参考程序
function experiment31 %观察单个正弦波泄漏、栅栏、正弦波的整周期取样
close all%关闭所有图形窗口
Ts=0.01;%取样周期
fs=1/Ts;
N=68;%取样点数68,取样长度为0.68秒，共6.8个信号周期
n=0:N-1;t=n*Ts;%计算取样时间
xn=cos(20*pi*t);%信号频率10Hz，周期0.1秒
subplot(221);stem(t,xn);%绘制原信号波形
axis([0,0.5,-10,10]);title(['T=0.01,N=',num2str(N)]);ylabel('x1(n)');
Xk=fft(xn);MAGXk=abs(Xk);%进行FFT变换并求幅频特性
fn=fs/N*n %对应各个n值的实际模拟频率值
subplot(222);stem(n,MAGXk);%绘制幅频特性
axis([0,20,0,50]);ylabel('X (k)');title('x(n)的频谱');
N=100;%取样点数100,取样长度为1秒，共10个信号周期
n=0:N-1;t=n*Ts;
xn=cos(20*pi*t);
subplot(223);stem(t,xn);
axis([0,0.5,-10,10]);title(['T=0.01,N=',num2str(N)]);ylabel('x2(n)');
Xk=fft(xn);MAGXk=abs(Xk);
fn=fs/N*n %对应各个n值的实际模拟频率值
subplot(224);stem(n,MAGXk);
ylabel('X(k)');title('x(n)的频谱');
function experiment32 %进一步观察泄漏及栅栏效应
close all%关闭所有图形窗口
Ts=0.01;%取样周期
fs=1/Ts;
N=100;%取样点数100,
n=0:N-1;t=n*Ts;%计算取样时间
xn=cos(20*pi*t)+2*cos(55*pi*t);%信号频率10Hz，27.5Hz,2倍幅值关系subplot(421);stem(t,xn);%绘制原信号波形
ylim([-3,3]);title(['T=0.005,N=',num2str(N)]);ylabel('x1(n)');
Xk=fft(xn);MAGXk=abs(Xk);%进行FFT变换并求幅频特性
fn=fs/N*n %对应各个n值的实际模拟频率值
subplot(422);stem(n,MAGXk);%绘制幅频特性
axis([0,50,0,100]);ylabel('X1(k)');title('x1(n)的频谱');
N=512;%取样点数扩大,取样长度增加
n=0:N-1;t=n*Ts;%计算取样时间
xn=cos(20*pi*t)+2*cos(55*pi*t);%信号频率10Hz，27.5Hz, 2倍幅值关系subplot(423);stem(t,xn);%绘制原信号波形
ylim([-3,3]);title(['T=0.005,N=',num2str(N)]);ylabel('x1(n)');
Xk=fft(xn);MAGXk=abs(Xk);%进行FFT变换并求幅频特性
fn=fs/N*n %对应各个n值的实际模拟频率值
subplot(424);stem(n,MAGXk);%绘制幅频特性
axis([0,256,0,512]);ylabel('X(k)');title('x(n)的频谱');
N=1024;%取样点数扩大,取样长度增加
n=0:N-1;t=n*Ts;%计算取样时间
xn=cos(20*pi*t)+2*cos(55*pi*t);
subplot(425);stem(t,xn);%绘制原信号波形
ylim([-3,3]);title(['T=0.005,N=',num2str(N)]);ylabel('x1(n)');
Xk=fft(xn);MAGXk=abs(Xk);%进行FFT变换并求幅频特性
fn=fs/N*n %对应各个n值的实际模拟频率值
subplot(426);stem(n,MAGXk);%绘制幅频特性
axis([0,512,0,1024]);ylabel('X(k)');title('x(n)的频谱');
N=2048;%取样点数扩大,取样长度增加
n=0:N-1;t=n*Ts;%计算取样时间
xn=cos(20*pi*t)+2*cos(55*pi*t);
subplot(427);stem(t,xn);%绘制原信号波形
ylim([-3,3]);title(['T=0.005,N=',num2str(N)]);ylabel('x1(n)');
Xk=fft(xn);MAGXk=abs(Xk);%进行FFT变换并求幅频特性
fn=fs/N*n %对应各个n值的实际模拟频率值
subplot(428);stem(n,MAGXk);%绘制幅频特性
axis([0,1024,0,2048]);ylabel('X(k)');title('x(n)的频谱');
function experiment33%观察谱间干扰；观察截断、栅栏效应的综合影响。
% 本实验事先已知频率。
% 但是：若事先不知频率，则N值需要取大一些，当改变N值时，若频谱变化不大，则可认为误差小到一定范围，截断、栅栏等影响减小到一定范围。
close all%关闭所有图形窗口
Ts=1;%取样周期1秒，频率1Hz
fs=1/Ts;
N=10;%取样点数,取样长度为NTs秒，频率分辨率为1/NTs Hz
n=0:N-1;t=n*Ts;%计算取样时间
xn=cos(0.48*pi*t)+cos(0.52*pi*t);%信号频率0.24Hz,0.26Hz
subplot(221);stem(t,xn);%绘制原信号波形
title(['T=1,N=',num2str(N)]);ylabel('x1(n)');
Xk=fft(xn);MAGXk=abs(Xk);%对x(n)进行FFT变换并求幅频特性
fn=fs/N*n %对应各个n值的实际模拟频率值
subplot(222);stem(n,MAGXk);%绘制x1(n)幅频特性
axis([0,N/2,0,max(MAGXk)]);
fk=n/N/Ts%观察相应的k值对应的实际模拟频率值
ylabel('X(k)');title('x(n)的频谱');
M=100;%该部分观察补零对于频谱的影响，频谱会更密，但无助消除频谱间的混淆。
m=0:M-1;
xn2=[xn,zeros(1,90)];
title(['T=1,N=',num2str(N)]);ylabel('x1(n)');
Xk2=fft(xn2);MAGXk2=abs(Xk2);%进行FFT变换并求幅频特性
subplot(223);stem(m,xn2);%绘制原信号波形
subplot(224);stem(m,MAGXk2);%绘制x1(n)幅频特性
axis([0,100/2,0,max(MAGXk2)]);
fm=n/M/Ts%观察相应的m值对应的实际模拟频率值
ylabel('X(k)');title('x(n)的频谱');
figure;
N=100;%可以清晰分辨两个频率，对于幅度的指示也对，由于N值的特殊性(24整周期和26整周期），避开了泄漏，FFT取值点在主峰和零点上。
n=0:N-1;t=n*Ts;%计算取样时间
xn=cos(0.48*pi*t)+cos(0.52*pi*t);%信号频率0.24Hz,0.26Hz
subplot(321);stem(t,xn);%绘制原信号波形
title(['T=1,N=',num2str(N)]);ylabel('x1(n)');
Xk=fft(xn);MAGXk=abs(Xk);%进行FFT变换并求幅频特性
subplot(322);stem(n,MAGXk);%绘制幅频特性
axis([0,N/2,0,max(MAGXk)+1]);
fn=fs/N*n %对应各个n值的实际模拟频率值
ylabel('X(k)');title('x(n)的频谱');
N=101;%取样点数,取样长度改变，幅度未能精确指示，两个频率泄漏后各自的频谱取样的位置对应的幅度不一样
n=0:N-1;t=n*Ts;%计算取样时间
xn=cos(0.48*pi*t)+cos(0.52*pi*t);%信号频率0.24Hz,0.26Hz
subplot(323);stem(t,xn);%绘制原信号波形
title(['T=1,N=',num2str(N)]);ylabel('x1(n)');
Xk=fft(xn);MAGXk=abs(Xk);%进行FFT变换并求幅频特性
subplot(324);stem(n,MAGXk);%绘制x1(n)幅频特性
axis([0,N/2,0,max(MAGXk)+1]);
fn=fs/N*n %对应各个n值的实际模拟频率值
ylabel('X(k)');title('x(n)的频谱');
N=102;%仍然不理想
n=0:N-1;t=n*Ts;%计算取样时间
xn=cos(0.48*pi*t)+cos(0.52*pi*t);%信号频率0.24Hz,0.26Hz
subplot(325);stem(t,xn);%绘制原信号波形
title(['T=1,N=',num2str(N)]);ylabel('x1(n)');
Xk=fft(xn);MAGXk=abs(Xk);%进行FFT变换并求幅频特性
subplot(326);stem(n,MAGXk);%绘制x1(n)幅频特性
axis([0,N/2,0,max(MAGXk)+1]);
fn=fs/N*n %对应各个n值的实际模拟频率值
ylabel('X(k)');title('x(n)的频谱');
figure;
N=512;%两频率幅度一致，周围泄漏有指示
n=0:N-1;t=n*Ts;%计算取样时间
xn=cos(0.48*pi*t)+cos(0.52*pi*t);%信号频率0.24Hz,0.26Hz
subplot(221);stem(t,xn);%绘制原信号波形
title(['T=1,N=',num2str(N)]);ylabel('x1(n)');
Xk=fft(xn);MAGXk=abs(Xk);%进行FFT变换并求幅频特性
subplot(222);stem(n,MAGXk);%绘制x1(n)幅频特性
axis([0,N/2,0,max(MAGXk)+1]);
fn=fs/N*n %对应各个n值的实际模拟频率值
ylabel('X(k)');title('x(n)的频谱');
N=513;%幅度不一致，周围泄漏有指示
n=0:N-1;t=n*Ts;%计算取样时间
xn=cos(0.48*pi*t)+cos(0.52*pi*t);%信号频率0.24Hz,0.26Hz
subplot(223);stem(t,xn);%绘制原信号波形
title(['T=1,N=',num2str(N)]);ylabel('x1(n)');
Xk=fft(xn);MAGXk=abs(Xk);%进行FFT变换并求幅频特性
subplot(224);stem(n,MAGXk);%绘制x1(n)幅频特性
实验二 IIR数字滤波器设计及软件实现 1．实验目的 （1）熟悉用双线性变换法设计IIR数字滤波器的原理与方法； （2）学会调用MATLAB信号处理工具箱中滤波器设计函数（或滤波器设计分析工具fdatool）设计各种IIR数字滤波器，学会根据滤波需求确定滤波器指标参数。 （3）掌握IIR数字滤波器的MATLAB实现方法。 （3）通过观察滤波器输入输出信号的时域波形及其频谱，建立数字滤波的概念。 2．实验原理 设计IIR数字滤波器一般采用间接法（脉冲响应不变法和双线性变换法），应用最广泛的是双线性变换法。基本设计过程是：①先将给定的数字滤波器的指标转换成过渡模拟滤波器的指标；②设计过渡模拟滤波器；③将过渡模拟滤波器系统函数转换成数字滤波器的系统函数。MATLAB信号处理工具箱中的各种IIR数字滤波器设计函数都是采用双线性变换法。第六章介绍的滤波器设计函数butter、cheby1 、cheby2 和ellip可以分别被调用来直接设计巴特沃斯、切比雪夫1、切比雪夫2和椭圆模拟和数字滤波器。本实验要求读者调用如上函数直接设计IIR数字滤波器。 本实验的数字滤波器的MATLAB实现是指调用MATLAB信号处理工具箱函数filter对给定的输入信号x(n)进行滤波，得到滤波后的输出信号y(n）。 3. 实验内容及步骤 （1）调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st，该函数还会自动绘图显示st的时域波形和幅频特性曲线，如图1所示。由图可见，三路信号时域混叠无法在时域分离。但频域是分离的，所以可以通过滤波的方法在频域分离，这就是本实验的目的。 图1 三路调幅信号st的时域波形和幅频特性曲线 （2）要求将st中三路调幅信号分离，通过观察st的幅频特性曲线，分别确定可以分离st中三路抑制载波单频调幅信号的三个滤波器（低通滤波器、带通滤波器、高通滤波器）的通带截止频率和阻带截止频率。要求滤波器的通带最大衰减为0.1dB,阻带最小衰减为
实验一MATLAB语言的基本使用方法 实验类别：基础性实验 实验目的： （1）了解MATLAB程序设计语言的基本方法，熟悉MATLAB软件运行环境。 （2）掌握创建、保存、打开m文件的方法，掌握设置文件路径的方法。 （3）掌握变量、函数等有关概念，具备初步的将一般数学问题转化为对应计算机模型并进行处理的能力。 （4）掌握二维平面图形的绘制方法，能够使用这些方法进行常用的数据可视化处理。 实验内容和步骤： 1、打开MATLAB，熟悉MATLAB环境。 2、在命令窗口中分别产生3*3全零矩阵，单位矩阵，全1矩阵。 3、学习m文件的建立、保存、打开、运行方法。 4、设有一模拟信号f(t)=1.5sin60πt,取?t=0.001,n=0,1,2,…,N-1进行抽样，得到 序列f(n),编写一个m文件sy1_1.m,分别用stem,plot,subplot等命令绘制32 点序列f(n)(N=32)的图形，给图形加入标注，图注，图例。 5、学习如何利用MATLAB帮助信息。 实验结果及分析： 1）全零矩阵 >> A=zeros(3,3) A = 0 0 0 0 0 0 0 0 0 2）单位矩阵 >> B=eye(3) B = 1 0 0 0 1 0 0 0 1 3）全1矩阵 >> C=ones(3) C = 1 1 1 1 1 1 1 1 1 4)sy1_1.m N=32; n=0:N-1; dt=0.001; t=n*dt; y=1.5*sin(60*pi*t); subplot(2,1,1), plot(t,y); xlabel('t'); ylabel('y=1.5*sin(60*pi*t)'); legend('正弦函数'); title('二维图形'); subplot(2,1,2), stem(t,y) xlabel('t'); ylabel('y=1.5*sin(60*pi*t)'); legend('序列函数'); title('条状图形'); 00.0050.010.0150.020.0250.030.035 t y = 1 . 5 * s i n ( 6 * p i * t ) 二维图形 00.0050.010.0150.020.0250.030.035 t y = 1 . 5 * s i n ( 6 * p i * t ) 条状图形
（1） 观察高斯序列的时域和幅频特性，固定信号)(n x a 中参数p=8,改变q 的 值，使q 分别等于2、4、8，观察他们的时域和幅频特性，了解当q 取不同值时，对信号序列的时域和幅频特性的影响；固定q=8，改变p,使p 分别等于8、13、14，观察参数p 变化对信号序列的时域和幅频特性的影响，注意p 等于多少时会发生明显的泄漏现象，混叠是否也随之出现？记录实验中观察到的现象，绘出相应的时域序列和幅频特性曲线。 ()() ?????≤≤=-其他0150,2n e n x q p n a 解：程序见附录程序一： P=8,q 变化时： t/T x a (n ) k X a (k ) t/T x a (n ) p=8 q=4 k X a (k ) p=8 q=4 t/T x a (n ) p=8 q=8 k X a (k ) p=8 q=8 幅频特性 时域特性
t/T x a (n ) p=8 q=8 k X a (k ) p=8 q=8 t/T x a (n ) 5 10 15 k X a (k ) p=13 q=8 t/T x a (n ) p=14 q=8 5 10 15 k X a (k ) p=14 q=8 时域特性幅频特性 分析： 由高斯序列表达式知n=p 为期对称轴； 当p 取固定值时，时域图都关于n=8对称截取长度为周期的整数倍，没有发生明显的泄漏现象；但存在混叠，当q 由2增加至8过程中，时域图形变化越来越平缓，中间包络越来越大，可能函数周期开始增加，频率降低，渐渐小于fs/2,混叠减弱； 当q 值固定不变，p 变化时，时域对称中轴右移，截取的时域长度渐渐地不再是周期的整数倍，开始无法代表一个周期，泄漏现象也来越明显，因而图形越来越偏离真实值， p=14时的泄漏现象最为明显，混叠可能也随之出现；
2.1 clc close all; n=0:15; p=8;q=2; x=exp(-(n-p.^2/q; figure(1; subplot(3,1,1; stem(n,x; title('exp(-(n-p^2/q,p=8,q=2'; xk1=fft(x,16; q=4; x=exp(-(n-p.^2/q; subplot(3,1,2; xk2=fft(x,16; stem(n,x; title('exp(-(n-p^2/q,p=8,q=4'; q=8; x=exp(-(n-p.^2/q;
xk3=fft(x,16; subplot(3,1,3; stem(n,x; title('exp(-(n-p^2/q,p=8,q=8';%时域特性figure(2; subplot(3,1,1; stem(n,abs(xk1; title('exp(-(n-p^2/q,p=8,q=2'; subplot(3,1,2; stem(n,abs(xk2; title('exp(-(n-p^2/q,p=8,q=4'; subplot(3,1,3; stem(n,abs(xk3; title('exp(-(n-p^2/q,p=8,q=8';%频域特性%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%% p=8;q=8; figure(3; subplot(3,1,1; stem(n,x; title('exp(-(n-p^2/q,p=8,q=8';
xk1=fft(x,16; p=13; x=exp(-(n-p.^2/q; subplot(3,1,2; xk2=fft(x,16; stem(n,x; title('exp(-(n-p^2/q,p=13,q=8'; p=14; x=exp(-(n-p.^2/q; xk3=fft(x,16; subplot(3,1,3; stem(n,x; title('exp(-(n-p^2/q,p=14,q=8';%时域特性figure(4; subplot(3,1,1; stem(n,abs(xk1; title('exp(-(n-p^2/q,p=8,q=8'; subplot(3,1,2; stem(n,abs(xk2; title('exp(-(n-p^2/q,p=13,q=8'; subplot(3,1,3;
实验6 数字滤波器的网络结构 一、实验目的： 1、加深对数字滤波器分类与结构的了解。 2、明确数字滤波器的基本结构及其相互间的转换方法。 3、掌握用MA TLAB 语言进行数字滤波器结构间相互转换的子函数及程序编写方法。 二、实验原理： 1、数字滤波器的分类 离散LSI 系统对信号的响应过程实际上就是对信号进行滤波的过程。因此，离散LSI 系统又称为数字滤波器。 数字滤波器从滤波功能上可以分为低通、高通、带通、带阻以及全通滤波器；根据单位脉冲响应的特性，又可以分为有限长单位脉冲响应滤波器（FIR ）和无限长单位脉冲响应滤波器（IIR ）。 一个离散LSI 系统可以用系统函数来表示： M -m -1-2-m m m=0 012m N -1-2-k -k 12k k k=1 b z b +b z +b z ++b z Y(z)b(z)H(z)=== =X(z)a(z) 1+a z +a z ++a z 1+a z ∑∑ 也可以用差分方程来表示： N M k m k=1 m=0 y(n)+a y(n-k)=b x(n-m)∑∑ 以上两个公式中，当a k 至少有一个不为0时，则在有限Z 平面上存在极点，表达的是以一个IIR 数字滤波器；当a k 全都为0时，系统不存在极点，表达的是一个FIR 数字滤波器。FIR 数字滤波器可以看成是IIR 数字滤波器的a k 全都为0时的一个特例。 IIR 数字滤波器的基本结构分为直接Ⅰ型、直接Ⅱ型、直接Ⅲ型、级联型和并联型。 FIR 数字滤波器的基本结构分为横截型（又称直接型或卷积型）、级联型、线性相位型及频率采样型等。本实验对线性相位型及频率采样型不做讨论，见实验10、12。 另外，滤波器的一种新型结构——格型结构也逐步投入应用，有全零点FIR 系统格型结构、全极点IIR 系统格型结构以及全零极点IIR 系统格型结构。 2、IIR 数字滤波器的基本结构与实现 （1）直接型与级联型、并联型的转换 例6-1 已知一个系统的传递函数为 -1-2-3 -1-2-3 8-4z +11z -2z H(z)=1-1.25z +0.75z -0.125z 将其从直接型（其信号流图如图6-1所示）转换为级联型和并联型。
实验一 MATLAB 仿真软件的基本操作命令和使用方法 实验容 1、帮助命令 使用 help 命令，查找 sqrt （开方）函数的使用方法； 2、MATLAB 命令窗口 （1）在MATLAB 命令窗口直接输入命令行计算3 1)5.0sin(21+=πy 的值； （2）求多项式 p(x) = x3 + 2x+ 4的根； 3、矩阵运算 （1）矩阵的乘法 已知 A=[1 2;3 4]， B=[5 5;7 8]，求 A^2*B
（2）矩阵的行列式 已知A=[1 2 3;4 5 6;7 8 9]，求A （3）矩阵的转置及共轭转置 已知A=[1 2 3;4 5 6;7 8 9]，求A' 已知B=[5+i,2-i,1;6*i,4,9-i], 求B.' , B' （4）特征值、特征向量、特征多项式 已知A=[1.2 3 5 0.9;5 1.7 5 6;3 9 0 1;1 2 3 4] ，求矩阵A的特征值、特征向量、特征多项式；
（5）使用冒号选出指定元素 已知：A=[1 2 3;4 5 6;7 8 9]；求A 中第3 列前2 个元素；A 中所有列第2，3 行的元素； 4、Matlab 基本编程方法 （1）编写命令文件：计算1+2+…+n<2000 时的最大n 值；
（2）编写函数文件：分别用for 和while 循环结构编写程序，求 2 的0 到15 次幂的和。
5、MATLAB基本绘图命令 （1）绘制余弦曲线 y=cos(t)，t∈[0，2π]
（2）在同一坐标系中绘制余弦曲线 y=cos(t-0.25)和正弦曲线 y=sin(t-0.5)， t∈[0，2π] （3）绘制[0，4π]区间上的 x1=10sint 曲线，并要求： （a）线形为点划线、颜色为红色、数据点标记为加号； （b）坐标轴控制：显示围、刻度线、比例、网络线 （c）标注控制：坐标轴名称、标题、相应文本； >> clear;
一、实验目的 1. 通过本次实验回忆并熟悉MATLAB这个软件。 2. 通过本次实验学会如何利用MATLAB进行序列的简单运算。 3. 通过本次实验深刻理解理论课上的数字信号处理的一个常见方法——对时刻n的样本附近的一些样本求平均，产生所需的输出信号。 3. 通过振幅调制信号的产生来理解载波信号与调制信号之间的关系。 二、实验内容 1. 编写程序在MATLAB中实现从被加性噪声污染的信号中移除噪声的算法，本次试验采用三点滑动平均算法，可直接输入程序P1.5。 2. 通过运行程序得出的结果回答习题Q1.31-Q1.33的问题，加深对算法思想的理解。 3. 编写程序在MATLAB中实现振幅调制信号产生的算法，可直接输入程序P1.6。 4. 通过运行程序得出的结果回答习题Q1.34-Q1.35的问题，加深对算法思想的理解。 三、主要算法与程序 1. 三点滑动平均算法的核心程序： %程序P1.5 %通过平均的信号平滑 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('时间序号n');ylabel('振幅'); legend('d[n]','s[n]','x[n]'); x1=[0 0 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('y[n]','s[n]'); xlabel('时间序号n');ylabel('振幅'); 2. 振幅调制信号的产生核心程序：（由于要几个结果，因此利用subplot函数画图） %程序P1.6 %振幅调制信号的产生 n=0:100; m=0.1;fH=0.1;fL=0.01; m1=0.3;fH1=0.3;fL1=0.03; xH=sin(2*pi*fH*n); xL=sin(2*pi*fL*n); y=(1+m*xL).*xH; xH1=sin(2*pi*fH1*n); xL1=sin(2*pi*fL1*n); y1=(1+m1*xL).*xH; y2=(1+m*xL).*xH1; y3=(1+m*xL1).*xH; subplot(2,2,1); stem(n,y); grid; xlabel('时间序号n');ylabel('振幅');title('m=0.1;fH=0.1;fL=0.01;'); subplot(2,2,2); stem(n,y1); grid; xlabel('时间序号n');ylabel('振幅');title('m=0.3;fH=0.1;fL=0.01;'); subplot(2,2,3); stem(n,y2); grid; xlabel('时间序号n');ylabel('振幅');title('m=0.3;fH=0.3;fL=0.01;'); subplot(2,2,4); stem(n,y3); grid;
武汉工程大学 数字信号处理实验报告 姓名：周权 学号：1204140228 班级：通信工程02
一、实验设备 计算机，MATLAB语言环境。 二、实验基础理论 1.序列的相关概念 2.常见序列 3.序列的基本运算 4.离散傅里叶变换的相关概念 5.Z变换的相关概念 三、实验内容与步骤 1.离散时间信号（序列）的产生 利用MATLAB语言编程产生和绘制单位样值信号、单位阶跃序列、指数序列、正弦序列及随机离散信号的波形表示。 四实验目的 认识常用的各种信号，理解其数字表达式和波形表示，掌握在计算机中生成及绘制数字信号波形的方法，掌握序列的简单运算及计算机实现与作用，理解离散时间傅里叶变换，Z变换及它们的性质和信号的频域分
实验一离散时间信号（序列）的产生 代码一 单位样值 x=2; y=1; stem(x,y); title('单位样值 ') 单位阶跃序列 n0=0; n1=-10; n2=10; n=[n1:n2]; x=[(n-n0)>=0]; stem(n,x); xlabel('n'); ylabel('x{n}'); title('单位阶跃序列');
实指数序列 n=[0:10]; x=(0.5).^n; stem(n,x); xlabel('n'); ylabel('x{n}'); title('实指数序列');
正弦序列 n=[-100:100]; x=2*sin(0.05*pi*n); stem(n,x); xlabel('n'); ylabel('x{n}'); title('正弦序列');
随机序列 n=[1:10]; x=rand(1,10); subplot(221); stem(n,x); xlabel('n'); ylabel('x{n}'); title('随机序列');
实验5 抽样定理 一、实验目的： 1、了解用MA TLAB 语言进行时域、频域抽样及信号重建的方法。 2、进一步加深对时域、频域抽样定理的基本原理的理解。 3、观察信号抽样与恢复的图形，掌握采样频率的确定方法和插公式的编程方法。 二、实验原理： 1、时域抽样与信号的重建 （1）对连续信号进行采样 例5-1 已知一个连续时间信号sin sin(),1Hz 3 ππ=0001f(t)=(2f t)+6f t f ，取最高有限带宽频率f m =5f 0，分别显示原连续时间信号波形和F s >2f m 、F s =2f m 、F s <2f m 三情况下抽样信号的波形。 程序清单如下： %分别取Fs=fm ，Fs=2fm ，Fs=3fm 来研究问题 dt=0.1; f0=1; T0=1/f0; m=5*f0; Tm=1/fm; t=-2:dt:2; f=sin(2*pi*f0*t)+1/3*sin(6*pi*f0*t); subplot(4,1,1); plot(t,f); axis([min(t),max(t),1.1*min(f),1.1*max(f)]); title('原连续信号和抽样信号'); for i=1:3; fs=i*fm;Ts=1/fs; n=-2:Ts:2; f=sin(2*pi*f0*n)+1/3*sin(6*pi*f0*n); subplot(4,1,i+1);stem(n,f,'filled'); axis([min(n),max(n),1.1*min(f),1.1*max(f)]); end 程序运行结果如图5-1所示：
原连续信号和抽样信号 图5-1 （2）连续信号和抽样信号的频谱 由理论分析可知，信号的频谱图可以很直观地反映出抽样信号能否恢复原模拟信号。因此，我们对上述三种情况下的时域信号求幅度谱，来进一步分析和验证时域抽样定理。 例5-2编程求解例5-1中连续信号及其三种抽样频率（F s>2f m、F s=2f m、F s<2f m）下的抽样信号的幅度谱。 程序清单如下： dt=0.1;f0=1;T0=1/f0;fm=5*f0;Tm=1/fm; t=-2:dt:2;N=length(t); f=sin(2*pi*f0*t)+1/3*sin(6*pi*f0*t); wm=2*pi*fm;k=0:N-1;w1=k*wm/N; F1=f*exp(-j*t'*w1)*dt;subplot(4,1,1);plot(w1/(2*pi),abs(F1)); axis([0,max(4*fm),1.1*min(abs(F1)),1.1*max(abs(F1))]); for i=1:3; if i<=2 c=0;else c=1;end fs=(i+c)*fm;Ts=1/fs; n=-2:Ts:2;N=length(n); f=sin(2*pi*f0*n)+1/3*sin(6*pi*f0*n); wm=2*pi*fs;k=0:N-1; w=k*wm/N;F=f*exp(-j*n'*w)*Ts; subplot(4,1,i+1);plot(w/(2*pi),abs(F)); axis([0,max(4*fm),1.1*min(abs(F)),1.1*max(abs(F))]); end 程序运行结果如图5-2所示。 由图可见，当满足F s≥2f m条件时，抽样信号的频谱没有混叠现象；当不满足F s≥2f m 条件时，抽样信号的频谱发生了混叠，即图5-2的第二行F s＜2f m的频谱图，，在f m=5f0的围，频谱出现了镜像对称的部分。
实验1认识实验-MATLAB语言上机操作实践 一、实验目的 ㈠了解MATLAB语言的主要特点、作用。 ㈡学会MATLAB主界面简单的操作使用方法。 ㈢学习简单的数组赋值、运算、绘图、流程控制编程。 二、实验原理 ㈠简单的数组赋值方法 MATLAB中的变量和常量都可以是数组（或矩阵）,且每个元素都可以是复数。 在MATLAB指令窗口输入数组A=[1 2 3；4 5 6；7 8 9]，观察输出结果。然后，键入：A（4，2）= 11 键入：A (5，：) = [-13 -14 -15] 键入：A（4，3）= abs (A(5，1)) 键入：A ([2，5]，：) = [ ] 键入：A/2 键入：A (4，：) = [sqrt(3) (4+5)/6*2 –7] 观察以上各输出结果。将A式中分号改为空格或逗号，情况又如何？请在每式的后面标注其含义。 2．在MATLAB指令窗口输入B=[1+2i，3+4i；5+6i ，7+8i], 观察输出结果。 键入：C=[1，3；5，7]+[2，4；6，8]*i，观察输出结果。 如果C式中i前的*号省略，结果如何？ 键入：D = sqrt (2+3i) 键入：D*D 键入：E = C’, F = conj(C), G = conj(C)’ 观察以上各输出结果, 请在每式的后面标注其含义。 3．在MATLAB指令窗口输入H1=ones(3,2)，H2=zeros(2,3)，H3=eye(4)，观察输出结果。 ㈡、数组的基本运算 1．输入A=[1 3 5]，B= [2 4 6]，求C=A+B，D=A-2，E=B-A 2．求F1=A*3，F2=A.*B，F3=A./B，F4=A.\B, F5=B.\A, F6=B.^A, F7=2./B, F8=B.\2 *3．求B'，Z1=A*B’，Z2=B’*A 观察以上各输出结果,比较各种运算的区别，理解其含义。 ㈢、常用函数及相应的信号波形显示 例1：显示曲线f(t)=2sin(2πt)，（t>0） ⅰ点击空白文档图标（New M-file），打开文本编辑器。 ⅱ键入：t=0:0.01:3; (1) f=2*sin(2*pi*t); (2) plot(t,f); title(‘f（t）-t曲线’)； xlabel（‘t’），ylabel（‘f（t）’）；
实验四 有限长单位脉冲响应滤波器设计 朱方方 0806020433 通信四班 (1) 设计一个线性相位FIR 高通滤波器，通带边界频率为0.6π，阻带边界频率为0.4π，阻 带衰减不小于40dB 。要求给出h(n)的解析式，并用MATLAB 绘出时域波形和幅频特性。 解： (1) 求数字边界频率: 0.6 , .c r ωπωπ== (2) 求理想滤波器的边界频率: 0.5n ωπ= (3) 求理想单位脉冲响应: []d s i n ()s i n [()] () ()1n n n n n n h n n παωαα παωα π?-- -≠??-=? ? -=?? (4) 选择窗函数。阻带最小衰减为-40dB ，因此选择海明窗(其阻带最小衰减为-44dB)；滤 波器的过渡带宽为0.6π-0.4π=0.2π，因此 6.21 0.231 , 152 N N N ππα-=?=== (5) 求FIR 滤波器的单位脉冲响应h(n): []31d sin (15)sin[0.5(15)] 1cos ()15()()()15(15)1 15 n n n R n n h n w n h n n n ππππ?---????-? ?≠? ???==-???? ? ?=? 程序： clear; N=31; n=0:N-1; hd=(sin(pi*(n-15))-sin(0.5*pi*(n-15)))./(pi *(n-15)); hd(16)=0.5; win=hanning(N); h=win'.*hd; figure; stem(n,h); xlabel('n'); ylabel('h(n)'); grid; title('FIR 高通滤波单位脉冲响应h(n)'); [H,w]=freqz(h,1); H=20*log10(abs(H)); figure;3 plot(w/pi,H); axis([0 1 -100 10]); xlabel('\omega/\pi'); ylabel('幅度/dB'); grid; title('FIR 高通滤波器，hanning 窗，N=31');
数字信号处理实验实验一离散时间信号与系统及MA TLAB实现 1.单位冲激信号： n = -5:5; x = (n==0); subplot(122); stem(n, x); 2.单位阶跃信号： x=zeros(1,11); n0=0; n1=-5; n2=5; n = n1:n2; x(:,n+6) = ((n-n0)>=0); stem(n,x); 3.正弦序列： n = 0:1/3200:1/100; x=3*sin(200*pi*n+1.2); stem(n,x); 4.指数序列 n = 0:1/2:10; x1= 3*(0.7.^n); x2=3*exp((0.7+j*314)*n); subplot(221); stem(n,x1); subplot(222); stem(n,x2); 5.信号延迟 n=0:20; Y1=sin(100*n); Y2=sin(100*(n-3)); subplot(221); stem(n,Y1); subplot(222); stem(n,Y2);
6.信号相加 X1=[2 0.5 0.9 1 0 0 0 0]; X2=[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7]; X=X1+X2; stem(X); 7.信号翻转 X1=[2 0.5 0.9 1]; n=1:4; X2=X1(5-n); subplot(221); stem(n,X1); subplot(222); stem(n,X2); 8.用MATLAB计算序列{-2 0 1 –1 3}和序列{1 2 0 -1}的离散卷积。a=[-2 0 1 -1 3]; b=[1 2 0 -1]; c=conv(a,b); M=length(c)-1; n=0:1:M; stem(n,c); xlabel('n'); ylabel('幅度'); 9．用MA TLAB计算差分方程 当输入序列为时的输出结果。 N=41; a=[0.8 -0.44 0.36 0.22]; b=[1 0.7 -0.45 -0.6]; x=[1 zeros(1,N-1)]; k=0:1:N-1; y=filter(a,b,x); stem(k,y) xlabel('n'); ylabel('幅度') 10.冲激响应impz N=64; a=[0.8 -0.44 0.36 0.22];
数字信处理上机实验答 案全 Document number【SA80SAB-SAA9SYT-SAATC-SA6UT-SA18】
第十章 上机实验 数字信号处理是一门理论和实际密切结合的课程，为深入掌握课程内容，最好在学习理论的同时，做习题和上机实验。上机实验不仅可以帮助读者深入的理解和消化基本理论，而且能锻炼初学者的独立解决问题的能力。本章在第二版的基础上编写了六个实验，前五个实验属基础理论实验，第六个属应用综合实验。 实验一系统响应及系统稳定性。 实验二时域采样与频域采样。 实验三用FFT对信号作频谱分析。 实验四 IIR数字滤波器设计及软件实现。 实验五 FIR数字滤波器设计与软件实现 实验六应用实验——数字信号处理在双音多频拨号系统中的应用 任课教师根据教学进度，安排学生上机进行实验。建议自学的读者在学习完第一章后作实验一；在学习完第三、四章后作实验二和实验三；实验四IIR数字滤波器设计及软件实现在。学习完第六章进行；实验五在学习完第七章后进行。实验六综合实验在学习完第七章或者再后些进行；实验六为综合实验，在学习完本课程后再进行。 实验一: 系统响应及系统稳定性 1.实验目的 （1）掌握求系统响应的方法。 （2）掌握时域离散系统的时域特性。 （3）分析、观察及检验系统的稳定性。 2.实验原理与方法 在时域中，描写系统特性的方法是差分方程和单位脉冲响应，在频域可以用系统函数描述系统特性。已知输入信号可以由差分方程、单位脉冲响应或系统函数求出系统对于该输入信号的响应，本实验仅在时域求解。在计算机上适合用递推法求差分方程的解，最简单的方法是采用MATLAB语言的工具箱函数filter函数。也可以用MATLAB语言的工具箱函数conv函数计算输入信号和系统的单位脉冲响应的线性卷积，求出系统的响应。 系统的时域特性指的是系统的线性时不变性质、因果性和稳定性。重点分析实验系统的稳定性，包括观察系统的暂态响应和稳定响应。 系统的稳定性是指对任意有界的输入信号，系统都能得到有界的系统响应。或者系统的单位脉冲响应满足绝对可和的条件。系统的稳定性由其差分方程的系数决定。 实际中检查系统是否稳定，不可能检查系统对所有有界的输入信号，输出是否都是有界输出，或者检查系统的单位脉冲响应满足绝对可和的条件。可行的方法是在系统的输入端加入单位阶跃序列，如果系统的输出趋近一个常数（包括零），就可以断定系统是稳定的[19]。系统的稳态输出是指当∞ n时，系统的输出。如果系统稳定，信号加入 → 系统后，系统输出的开始一段称为暂态效应，随n的加大，幅度趋于稳定，达到稳态输出。 注意在以下实验中均假设系统的初始状态为零。 3．实验内容及步骤
第二章上机作业 1、ljdt(A,B)函数定义 function ljdt(A,B) p=roots(A); q=roots(B); p=p'; q=q'; x=max(abs([p q 1])); x=x+0.1; y=x; clf hold on axis([-x x -y y]) w=0:pi/300:2*pi; t=exp(i*w); plot(t) axis('square') plot([-x x],[0 0]) plot([0 0],[-y y]) text(0.1,x,'jIm[z]') text(y,1/10,'Re[z]') plot(real(p),imag(p),'x') plot(ral(q),imag(q),'o') title('pole-zero diagram for discrete system') hold off 例2.26 a=[3 -1 0 0 0 1]; b=[1 1]; ljdt(a,b) p=roots(a) q=roots(b) pa=abs(p) 程序运行结果如下： P= 0.7255+0.4633i 0.7255+0.4633i -0.1861+0.7541i -0.1861-0.7541i -0.7455 q=
-1 pa= 0.8608 0.8608 0.7768 0.7768 0.7455 例2.27 b=[0 1 2 1];a=[1 -0.5 -0.005 0.3]; subplot 311 zplane(b,a);xlabel('实部');ylabel('虚部'); num=[0 1 2 1];den=[1 -0.5 -0.005 0.3]; h=impz(num,den); subplot 312
数字信号处理实验答案 HEN system office room 【HEN16H-HENS2AHENS8Q8-HENH1688】
实验一熟悉Matlab环境 一、实验目的 1.熟悉MATLAB的主要操作命令。 2.学会简单的矩阵输入和数据读写。 3.掌握简单的绘图命令。 4.用MATLAB编程并学会创建函数。 5.观察离散系统的频率响应。 二、实验内容 认真阅读本章附录，在MATLAB环境下重新做一遍附录中的例子，体会各条命令的含义。在熟悉了MATLAB基本命令的基础上，完成以下实验。 上机实验内容： （1）数组的加、减、乘、除和乘方运算。输入A=[1 2 3 4]，B=[3 4 5 6]，求C=A+B，D=A-B，E=A.*B，F=A./B，G=A.^B并用stem语句画出A、B、C、D、E、F、G。 clear all; a=[1 2 3 4]; b=[3 4 5 6]; c=a+b; d=a-b; e=a.*b; f=a./b; g=a.^b; n=1:4; subplot(4,2,1);stem(n,a); xlabel('n');xlim([0 5]);ylabel('A'); subplot(4,2,2);stem(n,b); xlabel('n');xlim([0 5]);ylabel('B'); subplot(4,2,3);stem(n,c); xlabel('n');xlim([0 5]);ylabel('C'); subplot(4,2,4);stem(n,d); xlabel('n');xlim([0 5]);ylabel('D'); subplot(4,2,5);stem(n,e); xlabel('n');xlim([0 5]);ylabel('E'); subplot(4,2,6);stem(n,f); xlabel('n');xlim([0 5]);ylabel('F'); subplot(4,2,7);stem(n,g); xlabel('n');xlim([0 5]);ylabel('G'); （2）用MATLAB实现下列序列： a) x(n)= 0≤n≤15 b) x(n)=e+3j)n 0≤n≤15 c) x(n)=3cosπn+π)+2sinπn+π) 0≤n≤15 d) 将c)中的x(n)扩展为以16为周期的函数x(n)=x(n+16),绘出四个周期。
文件名：tstem.m（实验一、二需要） 程序： f unction tstem(xn,yn) %时域序列绘图函数 %xn:被绘图的信号数据序列，yn:绘图信号的纵坐标名称（字符串）n=0:length(xn)-1; stem(n,xn,'.'); xlabel('n');ylabel('yn'); axis([0,n(end),min(xn),1.2*max(xn)]); 文件名：tplot.m（实验一、四需要） 程序： function tplot(xn,T,yn) %时域序列连续曲线绘图函数 %xn:信号数据序列,yn:绘图信号的纵坐标名称（字符串） %T为采样间隔 n=0;length(xn)-1;t=n*T; plot(t,xn); xlabel('t/s');ylabel(yn); axis([0,t(end),min(xn),1.2*max(xn)]); 文件名：myplot.m（实验一、四需要）
%(1)myplot;计算时域离散系统损耗函数并绘制曲线图。function myplot(B,A) %B为系统函数分子多项式系数向量 %A为系统函数分母多项式系数向量 [H,W]=freqz(B,A,1000) m=abs(H); plot(W/pi,20*log10(m/max(m)));grid on; xlabel('\omega/\pi');ylabel('幅度(dB)') axis([0,1,-80,5]);title('损耗函数曲线'); 文件名：mstem.m（实验一、三需要） 程序： function mstem(Xk) %mstem(Xk)绘制频域采样序列向量Xk的幅频特性图 M=length(Xk); k=0:M-1;wk=2*k/M;%产生M点DFT对应的采样点频率（关于pi归一化值） stem(wk,abs(Xk),'.');box on;%绘制M点DFT的幅频特性图xlabel('w/\pi');ylabel('幅度'); axis([0,2,0,1.2*max(abs(Xk))]); 文件名：mpplot.m（实验一需要）
数字信号处理作业实验题报告 第一章16.(1) 实验目的： 求解差分方程所描述的系统的单位脉冲响应和单位阶跃响应。 实验要求： 运用matlab求出y(n)=0.6y(n-1)-0.08y(n-2)+x(n)的单位脉冲响应和单位阶跃响应的示意图。 源程序： B1=1;A1=[1, -0.6, 0.08]; ys=2; %设差分方程 xn=[1, zeros(1, 20)]; %xn=单位脉冲序列，长度N=31 xi=filtic(B1, A1, ys); hn1=filter(B1, A1, xn, xi); %求系统输出信号hn1 n=0:length(hn1)-1; subplot(2, 1, 1);stem(n, hn1, '.') title('单位脉冲响应'); xlabel('n');ylabel('h(n)') xn=ones(1, 20); sn1=filter(B1, A1, xn, xi); %求系统输出信号sn1 n=0:length(sn1)-1; Subplot(2, 1, 2); stem(n, sn1, '.') title('单位阶跃响应'); xlabel('n'); ylabel('s(n)')
运行结果： 实验分析： 单位脉冲响应逐渐趋于0，阶跃响应保持不变，由此可见，是个稳定系统。
第二章31题 实验目的： 用matlab判断系统是否稳定。 实验要求： 用matlab画出系统的极，零点分布图，输入单位阶跃序列u(n)检查系统是否稳定。 源程序： A=[2, -2.98, 0.17, 2.3418, -1.5147]; B=[0, 0, 1, 5, -50]; subplot(2,1,1); zplane(B,A); %求H(z)的极点 p=roots(A); %求H(z)的模 pm=abs(p); if max(pm)<1 disp('系统因果稳定'), else,disp('系统因果不稳定'),end un=ones(1,800); sn=filter(B, A, un); n=0:length(sn)-1; subplot(2, 1, 2);plot(n, sn) xlabel('n');ylabel('s(n)')
§ Z 变换 ? Z 变换的定义及收敛域 【习题】 1. 假如)(n x 的z 变换代数表示式是下式，问)(z X 可能有多少不同的收敛域。 )83451)(411(411)(2122----+++- =z z z z z X 【分析】 ）要单独讨论，（环状、圆外、圆内：有三种收敛域：双边序列的收敛域为：特殊情况有：左边序列的收敛域为：因果序列的收敛域为：右边序列的收敛域为：特殊情况有：有限长序列的收敛域为 0 0 , , 0 0 , , 0 , 0 0 , 0 , 0 22 11 212 1∞==<<≤≤<≤<<≥≥∞≤<≥∞<<≤∞<≤≥∞≤<≤≤∞<<+ -++--z z R z R n n R z n n R z n n z R n n z R n z n z n n n z x x x x x x
解：对X (Z )的分子和分母进行因式分解得 )43 1 )(21 1)(211(2111111----+-+- =Z jZ jZ Z X (Z )的零点为：1/2，极点为：j/2，-j/2，-3/4 ∴ X (Z )的收敛域为： (1) 1/2 < | Z | < 3/4，为双边序列，见图一 (2) | Z | < 1/2，为左边序列，见图二 (3) | Z | > 3/4，为右边序列，见图三 图一 图二 图三 )431)(211)(411()211)(211()(11211-----++++- =Z Z Z Z Z Z X
? Z 反变换 【习题】 2. 有一右边序列 )(n x ，其 z 变换为)1)(211(1 )(11----=z z z X (a) 将上式作部分分式展开(用 1-z 表示)，由展开式求 )(n x 。 (b) 将上式表示成 z 的多项式之比，再作部分分式展开，由展开式求 )(n x ,并说明所得到的序列 与(a)所得的是一样的。 【注意】不管哪种表示法最后求出 x (n ) 应该是相同的。 解：(a) 因为11122 111)(---+--=z z z X 且x(n)是右边序列 所以 )()212()(n u n x n ?? ? ??-= (b) 122 1211 )1)(2 1(21231 )1)(2 1()(2 -+--+=---+=--=z z z z z z z z z X )()212( )1(2)1(21)()( n u n u n u n n x n n ??? ??-=-+-?? ? ??-=δ则
子程序： function myplot(B,A) %myplot(B,A) %时域离散系统损耗函数绘图 %B为系统函数分子多项式系数向量 %A为系统函数分母多项式系数向量 [H,W]=freqz(B,A,1000); m=abs(H); plot(W/pi,20*log10(m/max(m)));grid on; xlabel('\omega/\pi');ylabel('幅度(dB)') axis([0,1,-80,5]);title('损耗函数曲线'); function tplot(xn,T,yn) %时域序列连续曲线绘图函数 % xn:信号数据序列，yn:绘图信号的纵坐标名称（字符串） % T为采样间隔 n=0:length(xn)-1;t=n*T; plot(t,xn); xlabel('t/s');ylabel(yn); axis([0,t(end),min(xn),1.2*max(xn)]) 程序： %实验4程序exp4.m % IIR数字滤波器设计及软件实现 clear all;close all Fs=10000;T=1/Fs; %采样频率 %调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st st=mstg; %低通滤波器设计与实现========================================= fp=280;fs=450; wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60; %DF指标（低通滤波器的通、阻带边界频）[N,wp]=ellipord(wp,ws,rp,rs); %调用ellipord计算椭圆DF阶数N和通带截止频率wp [B,A]=ellip(N,rp,rs,wp); %调用ellip计算椭圆带通DF系统函数系数向量B和A y1t=filter(B,A,st); %滤波器软件实现 % 低通滤波器设计与实现绘图部分 figure(2);subplot(3,1,1); myplot(B,A); %调用绘图函数myplot绘制损耗函数曲线 yt='y_1(t)'; subplot(3,1,2);tplot(y1t,T,yt); %调用绘图函数tplot绘制滤波器输出波形 %带通滤波器设计与实现==================================================== fpl=440;fpu=560;fsl=275;fsu=900; wp=[2*fpl/Fs,2*fpu/Fs];ws=[2*fsl/Fs,2*fsu/Fs];rp=0.1;rs=60;