求两个序列的卷积和。
编写仿真程序，并调试得到结果，进行分析。
仿真程序及解释
nf1=0:
20;%定义步长为1的21点等间隔采样点，n为0～20
f1=0.8.^nf1;%序列f1的函数表达式f1=0.8n
subplot（2,2,1）;stem（nf1,f1,'filled'）;%视图以两行两列的形式显示，当前选定的窗体区域的序号为1，并定义横纵坐标及图表属性，画序列图
title（'f1（n）'）;%视图一标题为f1（n）
nf2=0:
10;%定义步长为1的11点等间隔采样点
lf2=length（nf2）;%求序列nf2的长度
f2=ones（1,lf2）;%序列nf2的函数表达式
subplot（2,2,2）;stem（nf2,f2,'filled'）;%视图以两行两列的形式显示，当前选定的窗体区域的序号为2，并定义横纵坐标及图表属性，画序列图
title（'f2（n）'）;%视图二标题为f2（n）
y=conv（f1,f2）;%求函数f1与f2的卷积和
subplot（2,1,2）;stem（y,'filled'）;title（'y（n）'）;%以两行一列的形式显示最后结果的图示，并定义横纵坐标及图表属性，视图标题为y（n）
调试结果及分析
✓线性时不变系统的输出等于输入序列和该系统的单位脉冲响应的卷积。
y（n）=∑0.8m*u（n-m）;卷积长度N=21+11-1=31，所以y（n）长度>=31.
实验二时域采样与频域采样
1.实验目的
（1）掌握模拟信号采样前后频谱的变化，以及如何选择采样频率才能使采样后的信号不丢失信息
（2）掌握频率域采样会引起时域周期化的概念，以及频率域采样定理及其对频域采样点数选择的指导作用。
（3）会用MATLAB语言进行时域抽样与信号重建的方法，以及频域抽样与恢复时程序的编写方法。
2.实验原理
了解时域采样定理的要点，理解理想采样信号
和模拟信号
之间的关系，
了解频域采样定理的要点，掌握这两个采样理论的结论：
“时域采样频谱周期延拓，频域采样时域信号周期延拓”。
3.实验内容
（1）时域采样理论的验证。
给定模拟信号，
式中A=444.128，
=50
π，
=50
πrad/s
（2）用DFT（FFT）求该模拟信号的幅频特性，选取三种采样频率，以验证时域采样理论。
（3）编写实验程序，计算
、
和
的幅度特性，并绘图显示。
观察分析频谱混叠失真。
仿真程序及解释
%时域采样理论验证
Tp=64/1000;
Fs=1000;T=1/Fs;
M=Tp*Fs;n=0:
M-1;%产生M=64采样序列x（n）
A=444.128;alph=pi*50*2^0.5;omega=pi*50*2^0.5;%给模拟信x（t）的系数赋值A=444.12,
=50
π，
=50
πrad/s
xnt=A*exp（-alph*n*T）.*sin（omega*n*T）;%进行nt采样，xnt=xa（nt）
Xk=T*fft（xnt,M）;%M点FFT[xnt],即64个点
yn=xnt;%赋值
subplot（3,2,1）;stem（n,yn,’k’）;boxon;title（'（a）Fs=1000Hz'）;%第一个图是xa（t）的Fs=1000Hz的时域抽样
k=0:
M-1;fk=k/Tp;
subplot（3,2,2）;plot（fk,abs（Xk））;title（'（a）T*FT[xa（nT）],Fs=1000Hz'）;%第二个图是xa（t）的幅频特性曲线
xlabel（'f（Hz）'）;ylabel（'幅度'）;axis（[0,Fs,0,1.2*max（abs（Xk））]）；%设置坐标比例
Fs=300;T=1/Fs;
M=ceil（Tp*Fs）;n=0:
M-1;%取整
xnt=A*exp（-alph*n*T）.*sin（omega*n*T）;%xnt=xa（nt）
Xk=T*fft（xnt,M）;%M点FFT[xnt）]
yn=xnt;
subplot（3,2,3）;stem（n,yn,’k’）;boxon;title（'（b）Fs=300Hz'）;%第三个图是xa（t）的Fs=300Hz的时域抽样
k=0:
M-1;fk=k/Tp;
subplot（3,2,4）;plot（fk,abs（Xk））;title（'（b）T*FT[xa（nT）,Fs=300Hz'）;%第四个图是xa（t）的幅频特性曲线
xlabel（'f（Hz）'）;ylabel（'幅度'）;axis（[0,Fs,0,1.2*max（abs（Xk））]）；%设置坐标比例
Fs=200;T=1/Fs;
M=ceil（Tp*Fs）;n=0:
M-1;
xnt=A*exp（-alph*n*T）.*sin（omega*n*T）;%xnt=xa（nt）
Xk=T*fft（xnt,M）;%M点FFT[xnt]
yn=xnt;subplot（3,2,5）;stem（n,yn,’k’）;boxon;title（'（c）Fs=200Hz'）;%第五个图是xa（t）的Fs=200Hz的时域抽样
k=0:
M-1;fk=k/Tp;
subplot（3,2,6）;plot（fk,abs（Xk））;title（'（c）T*FT[xa（nT）,Fs=200Hz'）;%第六个图是xa（t）的幅频特性曲线
xlabel（'f（Hz）'）;ylabel（'幅度'）;axis（[0,Fs,0,1.2*max（abs（Xk））]）；%设置坐标比例
调试结果及分析
式中A=444.128，
=50
π，
=50
πrad/s
✓^xa（t）=∑xa（t）δ（t-nT）,pδ（t）=∑δ（t-nT），其中n=-∞~∞，δ（t）是单位冲激信号，只有当t=nT时，才有非零值，故^xa（t）=∑xa（nT）δ（t-nT）
在傅里叶变换中，两信号在时域相乘的傅里叶变换等于两个信号分别的傅里叶变换的卷积。
ˆXa（jΩ）=Xa（jΩ）*Pδ（jΩ）／2π==∑xa（nT）e^（-jΩnt）
在数值上xa（nT）=x（n）故ˆXa（jΩ）=∑x（n）e^（-jωn），n=-∞~∞，ω=ΩT
✓理想采样信号的频谱是原模拟信号的频谱以Ωs为周期，进行周期性延拓而产生的；采样频率Ωs必须大于等于模拟信号最高频率两倍以上，才能使采样信号的频谱不产生频谱混叠。
在采样频率不变的情况下，采样点数越多，频谱分辨率越高。
如图当采样频率为1000Hz时，频谱混叠很小；当采样频率为300Hz时，频谱混叠严重；当采样频率为200Hz时，频谱混叠很严重。
如图所示，频率在0～500，0～150，0～100时周期再现。
（4）频域采样理论的验证
给定信号如下：
（5）编写程序分别对频谱函数
在区间
上等间隔采样32
和16点，得到
，再分别对
进行32点和16点IFFT，得到
。
仿真程序及解释
%频域采样理论验证
M=27;N=32;n=0:
M;
%产生M长三角波序列x（n）
xa=0:
floor（M/2）;xb=ceil（M/2）-1:
-1:
0;xn=[xa,xb];
Xk=fft（xn,1024）;%1024点FFT[x（n）],用于近似序列x（n）的TF
X32k=fft（xn,32）;%32点FFT[x（n）]
x32n=ifft（X32k）;%32点IFFT[X32（k）]得到x32（n）
X16k=X32k（1:
2:
N）;%隔点抽取X32k得到X16（K）
x16n=ifft（X16k,N/2）;%16点IFFT[X16（k）]得到x16（n）
subplot（3,2,2）;stem（n,xn,'.'）;boxon
title（'（b）三角波序列x（n）'）;
xlabel（'n'）;ylabel（'x（n）'）;axis（[0,32,0,20]）；%横纵轴的长度
k=0:
1023;%取样点0~1023
wk=2*k/1024;%赋值
subplot（3,2,1）;plot（wk,abs（Xk））;title（'（a）FT[x（n）]'）;%视图显示为三行两列，当前选定的窗体区域序号为1
xlabel（'\omega/\pi'）;ylabel（'|X（e^j^\omega）|'）;
%x轴y轴内容
axis（[0,1,0,200]）；%横纵轴的取值范围
k=0:
N/2-1;
subplot（3,2,3）;stem（k,abs（X16k）,'.'）;boxon
title（'（c）16点频域采样'）;
xlabel（'k'）;ylabel（'|X_1_6（k）|'）;axis（[0,8,0,200]）
n1=0:
N/2-1;
subplot（3,2,4）;stem（n1,x16n,'.'）;boxon
title（'（d）16点IDFT[X_1_6（k）]'）;xlabel（'n'）;ylabel（'x_1_6（n）'）;axis（[0,32,0,20]）
k=0:
N-1;
subplot（3,2,5）;stem（k,abs（X32k）,'.'）;boxon
title（'（e）32点频域采样'）;xlabel（'k'）;ylabel（'|X_3_2（k）|'）;axis（[0,16,0,200]）
n1=0:
N-1;
subplot（3,2,6）;stem（n1,x32n,'.'）;boxon
title（'（f）32点IDFT[X_3_2（k）]'）;xlabel（'n'）;ylabel（'x_3_2（n）'）;axis（[0,32,0,20]）
调试结果及分析
✓频域采样时域信号周期延拓，由于实序列的DFT满足共轭对称性，因此频域图仅画在[0,π]上的幅频特性波形。
✓频率采样点数N必须大于等于时域离散信号的长度M（即N≥M），才能使时域不产生混叠，则N点IDFT[XN（k）]得到的序列xN（n）就是原序列x（n）,即xN（n）=x（n）。
如果N>M,xN（n）比原序列尾部多N-M个零点；如果N如图所示，当采样点数N=16M时，无时域混叠失真，x32（n）=IDFT[X32（k）]=x（n）。
（6）分别画出
、
的幅度谱，并绘图显示x（n）、
的波形，进行对比和分析，验证总结频域采样理论。
4.思考题：
如果序列x（n）的长度为M，希望得到其频谱
在
上的N点等间隔采样，当N答：
先对原序列x（n）以N为周期进行周期延拓后取主值区序列，x（n）=[∑x（n+iN）]R（n）,再计算N点DFT，则得到N点频域采样：
XN（k）=DFT[x（n）]=
|w=2πk/Nk=0,1…N-1
实验三用双线性变换法设计IIR数字滤波器
1.实验目的
（1）熟悉用双线性变换法设计IIR数字滤波器的原理与方法。
（2）掌握数字滤波器的计算机仿真方法。
（3）通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。
2.实验内容
（1）用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。
设计指标参数为：
在通带内频率低于0.2π时，最大衰减小于1dB；在阻带内[0.3π,π]频率区间上，最小衰减大于15dB。
（2）以0.02π为采样间隔，打印出数字滤波器在频率区间［0,π/2］上的幅频响应特性曲线。
（3）用所设计的滤波器对实际心电图信号采样序列（在本实验后面给出）进行仿真滤波处理，并分别打印出滤波前后的心电图信号波形图，观察总结滤波作用与效果。
3.实验步骤
（1）复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR数字滤波器的内容，按照例6.4.2，用双线性变换法设计数字滤波器系统函数H（z）。
例6.4.2中已求出满足本实验要求的数字滤波器系统函数：
（2－1）
（2－2）
A=0.09036
B1=1.2686,C1=－0.7051
B2=1.0106,C2=－0.3583
B3=0.9044,C3=－0.2155
由（2－1）式和（2－2）式可见，滤波器H（z）由三个二阶滤波器H1（z），H2（z）和H3（z）级联组成，如图2－1所示。
图2－1滤波器H（z）的组成
（2）编写滤波器仿真程序,计算H（z）对心电图信号采样序列x（n）的响应序列y（n）。
设为第k级二阶滤波器的输出序列，为输入序列，如图2-1所示。
由（2-2）式可得到差分方程：
（2-3）
当k=1时，所以的总响应序列
可以用顺序迭代算法得到。
即依次对k=1,2,3,求解差分方程（2-3）,
最后得到仿真程序就是实现上述求解差分方程和顺序迭代算法的通用程序。
也可以直接调用MATLABfilter函数实现仿真。
（3）在通用计算机上运行仿真滤波程序，并调用通用绘图子程序，完成实验内容
（2）和（3）。
仿真程序及解释
x=[-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];
k=1;%给k赋初值，使滤波从第一级开始
closeall;%关闭所有参数
figure
（1）;%图标
subplot（2,2,1）;%将x（n）序列用图形显示出来
n=0:
55;%n值从0取到55
stem（n,x,'.'）;%画出序列函数离散点为n，幅度响应为x
axis（[055-10050]）;%纵坐标从0到50；横坐标从-100到50
holdon;%保存当前图像不被覆盖
n=0:
60;%n值从0取到60
m=zeros（61）;%提供零基准线
plot（n,m）;%画平面图
xlabel（'n'）;%x坐标显示n
ylabel（'x（n）'）;%y坐标显示x（n）
title（'心电图信号采样序列x（n）'）;
B=[0.090362*0.090360.09036];%输入序列的增益系数
A=[1.2686-0.7051];%第一级滤波器输出序列的增益系数
A1=[1.0106-0.3583];%第二级滤波器输出序列的增益系数
A2=[0.9044-0.2155];%第三级滤波器输出序列的增益系数
while（k<4）
y=filter（B,A,x）;
x=y;
if（k==2）
A=A1;
elseif（k==3）
A=A2;
else;
end
k=k+1;
end;%完成三级滤波
subplot（2,2,3）;%绘图
n=0:
55;%n值从0取到55
stem（n,y,'.'）;%画出序列函数离散点为n，幅度响应为y
axis（[055-155]）;%纵坐标从0到55；横坐标从-15到5
holdon;%保存当前图像不被覆盖
n=0:
60;%n值从0取到60
m=zeros（61）;%提供零基准线
plot（n,m）;%画平面图
xlabel（'n'）;%x坐标显示n
ylabel（'y（n）'）;%y坐标显示y（n）
title（'三级滤波后的心电图信号'）;%显示“三级滤波后的心电图信号”的图题
A=[0.09036,0.1872,0.09036];%输入序列的增益系数
B1=[1,-1.2686,0.7051];%第一级滤波器输出序列的增益系数
B2=[1,-1.0106,0.3583];%第二级滤波器输出序列的增益系数
B3=[1,-0.9044,0.2155];%第三级滤波器输出序列的增益系数
[H1,w]=freqz（A,B1,100）;%求第一级滤波器的频率响应
[H2,w]=freqz（A,B2,100）;%求第二级滤波器的频率响应
[H3,w]=freqz（A,B3,100）;%求第三级滤波器的频率响应
H4=H1.*（H2）;%第三级的频率响应等于前两级频率响应的卷积
H=H4.*（H3）;%一级二级滤波器矩阵相乘结果与三级滤波器矩阵相乘
mag=abs（H）;%求幅度变化
db=20*log10（（mag+eps）/max（mag））;%通带最大衰减幅度
subplot（2,2,2）%绘图
plot（w/pi,db）;%画平面图
axis（[00.5-5010]）;%纵坐标从0到0.5；横坐标从-50到10
title（'滤波器的幅频响应曲线'）;%显示“滤波器的幅频响应曲线”的图题
运行结果及分析
✓如图所示，图中第1和第3张图中红色部分是由程序中n=0:
60和m=zeros（61）两句实现；
图中第2张图中0.2表示通带截止频率；
程序中的while循环部分完成三级滤波的作用；
程序中A=[0.09036,0.1872,0.09036]中的0.09036是式子开三次方所得。
✓双线性变换法的特点：
采用非线性频率压缩的方法，将整个模拟频率轴压缩到±π/T之间，再变换到Z平面内，使之不产生频谱混叠现象，但是也不能保真的模仿模拟滤波器的频响曲线形状。
比较图中第1和第3张图
升级会员 