实验一 信号、系统及系统响应
一、实验目的
1、 熟悉连续信号经理想采样前后的频谱变化关系,加深对时域采样定理的理解。
2、 熟悉时域离散系统的时域特性。
3、 利用卷积方法观察分析系统的时域特性。
4、 掌握序列傅里叶变换的计算机实现方法,利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。
二、实验内容及步骤
1、 认真复习采样理论、离散信号与系统、线性卷积、序列的傅里叶变换及性质等有关内容,阅读本实验原理与方法。
2、 了解要使用到的MATLAB命令:
1)基于DTFT离散时间信号分析函数:freqz,real,imag,abs,angle,unwrap。
函数freqz可以用来计算一个以的有理分式形式给出的序列的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<N,则x先进行补零扩展为N点序列再求N点DFT。函数ifft(X)可以计算R点的谱序列的R点IDFT值,而ifft(X,N)同fft(x,N)的情况。
3)卷积函数conv用于求两个有限长度序列的卷积。
3、 实验准备工作:
1)产生实验中要用到的下列信号序列:
a.采样信号序列:对下面连续信号
进行采样,可得到采样序列
,
其中为幅度因子,为衰减因子,是模拟角频率,为采样间隔。这些参数都要在实验过程中由键盘输入,产生不同的和。
b.单位脉冲序列:
c.矩形序列: ,
2)产生系统单位脉冲响应序列。
a.
b.
3)编写有限长序列线性卷积子程序,用于完成两个给定长度的序列的卷积。也可以直接调用MATLAB语言中的卷积函数conv。conv用于两个有限长度序列的卷积,它假定两个序列都从n=0开始。调用格式如下:
y=conv(x,h)
其中参数x和h是两个已赋值的行向量序列。
4、完成下述实验内容:
1)分析采样序列的特性。产生采样序列,,,。
a、 取采样频率,即。观察所采样的幅频特性和的幅频特性在折叠频率处有无明显差别。应当注意,实验中所得频谱是用序列的傅立叶变换公式求得的,所以在频率量度上存在关系:。
b、 改变采样频率,,观察的变化并做记录。
c、 进一步降低采样频率,,观察频谱混叠是否明显存在,说明原因,并记录。
2) 离散信号、系统和系统响应分析。
a、观察信号和系统的时域和频域持性;利用线性卷积求信号通过系统的响应y(n),比较所求响应y(n)和的时域及频域特性,注意它们之间有无差异,绘图说明,并用所学结论解释所得结果。
b、 观察系统对信号的响应特性。利用线性卷积求系统响应y(n),并判断y(n)图形及其非零值序列长度是否与理论结果一致,对,说出一种定性判断y(n)图形是否正确的方法。调用序列傅立叶变换子程序,求得,观察的特性曲线,定性判断结果的正确性。改变的长度,取N=5,重复该试验。注意参数变化的影响,说明变化前后的差异,并解释所得的结果。
3) 卷积定理验证。将实验2)中的信号换成,使a=0.4,,A=1,T=1,重复实验2)a,打印曲线;对主程序作简单修改,计算,并绘出曲线,与前面直接对进行傅立叶变换所得幅频特性曲线进行比较,验证时域卷积定理。
三、思考题
1、 在分析理想采样序列特性的实验中;采样频率不同时,相应理想采样序列的傅里叶变换谱的数学频率度量是否都相同?它们所对应的模拟频率是否相同?为什么?
2、 在卷积定理验证的实验中,如果选用不同的频域采样点数M值,例如,选M=10和M=20,分别做序列的傅里叶变换,求得
,=0,1,……M-1
所得结果之间有无差异?为什么?
四、实验报告要求
1、 简述实验目的及实验原理。.
2、 按实验步骤附上实验过程中的信号序列、系统单位脉冲响应及系统响应序列的时域和幅频特性曲线,并对所得结果进行分析和解释。
3、 总结实验中的主要结论。
4、 简要回答思考题。
实验二 用FFT作谱分析
一、实验目的
1、 进一步加深DFT算法原理和基本性质的理解(因为FFT只是DFT的一种快速算法,所以FFT的运算结果必然满足DFT的基本性质)。
2、 熟悉FFT算法原理和FFT子程序的应用。
3、 学习用FFT对连续信号和时域离散信号进行谱分析的方法,了解可能出现的分析误差及其原因,以便在实际中正确应用FFT。
二、实验内容及步骤
1、 复习DFT的定义、性质和用DFT作谱分析的有关内容。
2、 复习FFT算法原理与编程思想,并对照DIT—FFT运算流图和程序框图,读懂本实验提供的FFT子程序。
3、 编制信号产生子程序,产生以下典型信号供谱分析用:
应当注意,如果给出的是连续信号,则首先要根据其最高频率确定采样频率,以及由频率分辨力确定采样点数N,然后对其进行软件采样(即计算,),产生对应序列。对信号,频率分辨力的选择要以能分辨开三个频率对应的谱线为准则。对周期序列,最好截取周期的整数倍进行谱分析,否则可能产生较大的分析误差,请同学们根据DFT的隐含周期性思考这个问题。
4、完成下述实验内容:
1) 对2中所给出的信号逐个进行谱分析。下面给出各个信号的FFT变换区间N以及连续信号的采样频率,供实验时参考。
:N=8,64
,:N=8,16
,:N=8,16
2)令,用FFT计算8点和16点离散傅里叶变换。
3)令,用FFT计算8点和16点离散傅里叶变换。
三、思考题
1、在N=8时,和的幅频特性会相同吗?为什么?N=16呢?
2、如果周期信号的周期预先不知道,如何使用FFT进行谱分析?
四、实验报告要求
1、简述实验原理及目的。
2、结合实验中所得给定典型序列幅频特性曲线,与理论结果比较,并分析说明误差产生的原因以及用FFT作谱分析时有关参数的选择方法。
3、总结实验所得主要结论。
4、简要回答思考题。
实验三 用双线性变换法设计IIR数字滤波器
一、实验目的
1、熟悉用双线性变换法设计IIR数字滤波器的原理与方法。
2、掌握数字滤波器的计算机仿真方法。
3、通过观察对实际心电图信号的滤波作用,获得数字滤波的感性知识。
二、实验内容
1、 用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。设计指标参数为:在通带内频率低于0.2时,最大衰减小于1dB;在阻带内[0.3,]频率区间上,最小衰减大于15dB。
2、 以T=1s为采样间隔,打印出数字滤波器在频率区间[0,]上的幅频响应特性曲线。
3、 用所设计的滤波器对实际心电图信号采样序列(在本实验后面给出)进行仿真滤波处理,并分别打印出滤波前后的心电图信号波形图,观察总结滤波作用与效果。
三、实验步骤
1、 复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR数字滤波器的内容,用双线性变换法设计数字滤波器系统函数。
教材中满足本实验要求的数字滤波器函数
式中:,k=1,2,3
A=0.09036,
B1=1.2686, C1=-0.7051
B2=1.0106, C2=-0.3583
B3=0.9044, C3=-0.2155
根据设计指标,调用MATLAB信号处理工具函数buttord和butter,也可得到。
由滤波器的函数可见,滤波器由三个二阶滤波器、、级联组成。
2、 编写滤波器仿真程序,计算对心电图信号采样序列 的响应序列y(n)。
为第k阶滤波器的输出序列,为输入序列,根据滤波器的组成可得差分方程
当时,。所以对的总响序列可以用顺序迭代算法得到。即依次对k=1,2,3求解差分方程,最后得到。仿真程序就是实现上述求解差分方程和顺序迭代算法的通用程序,也可直接调用MATLAB filter函数实现仿真。
3、 在通用计算机上运行仿真滤波程序,完成实验内容(2)和(3)。
四、实验报告要求
1、 简述实验目的及原理。
2、 由所打印的特性曲线及设计过程简述双线性变换法的特点。
五、心电图信号采样序列:
人体心电图信号在测量过程中往往受到工业高频干扰,所以必须经过低通滤波处理后,才能作为判断心脏功能的有用信息。下面给出一实际心电图信号采样序列样本,其中存在高频干扰。在实验中,以作为输入序列,滤除其中的干扰成分。
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}
实验四 用窗函数法设计FIR数字滤波器
一、实验目的:
1、掌握用窗函数法设计FIR数字滤波器的原理和方法。
2、熟悉线性相位FIR数字滤波器特性。
3、了解各种窗函数对滤波特性的影响。
二、实验内容及步骤:
1、复习用窗函数法设计FIR数字滤波器一节内容。
一般,设计线性相位FIR数字滤波器采用窗函数法或频率抽样法,本实验采用窗函数法,采用矩形窗。下面简要介绍窗函数法设计FIR理论知识,更详细的相关知识请参考教材。
如果所希望的滤波器的理想频率响应函数为,如理想的低通,由信号系统的知识知道,在时域系统的冲激响应hd(n)将是无限长的,如图1、图2所示。
图1 图2
若时域响应是无限长的,则不可能实现,因此需要对其截断,即设计一个FIR滤波器频率响应来逼近,即用一个窗函数w(n)来截断hd(n),如式3所示:
(式3)。
最简单的截断方法是矩形窗,实际操作中,直接取hd(n)的主要数据即可。
作为实际设计的FIR数字滤波器的单位脉冲响应序列,其频率响应函数为:
(式4)
令,则
(式5),式中,N为所选窗函数的长度。
如果要求线性相位特性,还必须满足:
(式6),根据式6中的正、负和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波器特性正确选择其中一类。例如:要设计线性相位低通特性,可选择类。
有关窗函数的种类及特性,请参照教材。
2、编写程序
1)编写能产生矩形窗、升余弦窗、改进升余弦窗和二阶升余弦窗的窗函数子程序,也可调用MATLAB相关函数。
2)编写低通FIR数字滤波器主程序。
3)上机实验内容。
A、用升余弦窗设计一线性相位低通FIR数字滤波器,截止频率。窗口长度N=15,33。要求在两种窗口长度情况下,分别求出(n),打印出相应的幅频特性和相频特性曲线,观察3dB带宽和20dB带宽。总结窗口长度N对滤波特性的影响。
B、=33,=,用四种窗函数设计线性相位低通滤波器。绘制相应的幅频特性曲线,观察3dB和20dB带宽以及阻带最小衰减,比较四种窗函数对滤波器特性的影响。
三、思考题:
1、 如果给定通带截止频率和阻带截止频率以及阻带最小衰减,如何用窗函数法设计线性相位低通滤波器?写出设计步骤。
2、 如果要求用窗函数法设计带通滤波器,且给定上、下边带截止频率为和,试求理想带通的单位脉冲响应。
四、实验报告要求:
1、 简述实验目的及原理。
2、 按照实验步骤及要求,比较各种情况下的滤波性能,说明窗口长度N和窗函数类型对滤波特性的影响。
3、 总结用窗函数法设计FIR滤波器的主要特点。
4、 简要回答思考题。
参考:用于数字信号处理的MATLAB操作或者函数
1、向量的生成:利用冒号(:)生成向量。
1)X=j:k (j<k)
生成X=[j,j+1,j+2,…,k-1,k]
2)X=j:I:k
如果I>0 且j<k或I<0且j>k,则生成向量
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 是同维矩阵对应位置元素相除;
表示矩阵的转置运算;
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
本实验可能用到的一些序列和函数
1)单位冲击序列:
function [x,n]=impseq(n0,n1,n2)
n=[n1:n2];
x=[(n-n0)==0];利用逻辑运算实现在n=N处为1。如:impseq(5,0,10)
2)单位阶跃序列:
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)实现了
也可以用
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*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)
函数调用
常用的函数: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)
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