数字信号处理实验报告
数字信号处理实验报告 姓 名:
班 级:
13 电信 2
学 号:
2013302
2013302
指导老师:
日期:
2016、6、6～2016、6、17
华南农业大学电子工程学院 电子信息工程系 实验一
常 见 离散信号的 B MATLAB 产生与图形显示
一、实验目的
加深对常用离散信号的理解; 二、实验原理
1． 单位抽样序列 01) ( n
00nn
在 MATLAB 中可以利用 zeros()函数实现。
; 1 ) 1 (); , 1 (xN zeros x 如果 ) (n 在时间轴上延迟了 k 个单位,得到 ) ( k n 即: 01) ( k n
0 nk n
2.
单位阶越序列
01) (n u
00nn
在 MATLAB 中可以利用 ones()函数实现。
); , 1 ( N ones x
3.
正弦序列
) / 2 sin( ) ( Fs fn A n x
在 MATLAB 中
) / * * * 2 sin( *1 : 0fai Fs n f pi A xN n 4.
复正弦序列
n je n x ) (
在 MATLAB 中
) * * exp(1 : 0n w j xN n 5.
指数序列
na n x ) (
在 MATLAB 中
n a xN n.^1 : 0 6、卷积分析 conv、m 用来实现两个离散序列的线性卷积。
其调用格式就是:y=conv(x,h) 若 x 的长度为 N,h 的长度为 M,则 y 的长度 L=N+M-1。
三、实验内容
1、画出信号 x(n) = 1、5*(n+1) - (n-3)的波形。
2、求序列 x(n)与 h(n)的线性卷积 y(n)=x(n)*h(n)。x(n) = {3,-3,7,0,-1,5,2} , h(n) = {2,3,0,-5,2,1}、 画出 x(n),h(n),y(n)与 n 的离散序列图形 四、实验要求 1)画出信号 x(n) = 1、5*(n+1) - (n-3)的波形。
①B MATLAB 程序 如下: n3 = [-3:3]; x3 = [(n3+1)==0];
subplot(1,3,1); stem(n3,x3); n4 = [-3:3]; x4 = [(n4-3)==0]; subplot(1,3,2);stem(n4,x4); n5 = [-3:3];x5 = 1、5*x3 - x4; subplot(1,3,3);stem(n5,x5); ②理论 计算: x(n)=
③程序 运行结果: 图(1) 从图(1)左侧起第一幅图就是信号(n+1)的波形,第二幅图就是信号(n-3)的波形,最后一幅图就是信号 x(n) = 1、5*(n+1) - (n-3)的波形。
2)求序列 x(n)与 h(n)的线性卷积 y(n)=x(n)*h(n);x(n) = {3,-3,7,0,-1,5,2} , h(n) = {2,3,0,-5,2,1},画出 x(n),h(n),y(n)与 n 的离散序列图形。
①B MATLAB 程序 如下: n6 = [0:6];x6 = [3,-3,7,0,-1,5,2];subplot(1,3,1);stem(n6,x6); n7 = [0:5];x7 = [2,3,0,-5,2,1];subplot(1,3,2);stem(n7,x7);n8 =
[0:11]; x8 = conv(x6,x7);subplot(1,3,3);stem(n8,x8); ②理论 分析: 信号 x(n)的长度为 7,,h(n)的长度为 6,则线性卷积 y(n)=x(n)*h(n)的长度为 7+6-1=12。
y(n) = {6 ,3 ,5 ,6,-7 ,-25 ,30 ,21 ,-23 ,-1 ,9 ,2} ③程序 运行结果: 图(2) 从图(2)左侧起第一幅图就是信号 x(n)的波形,第二幅图就是信号 h(n)的波形,最后一幅图就是线性卷积信号 y(n)=x(n)*h(n)的波形。经过比较,理论与实验结果一致。
实验二
离散系统的差分方程、冲激响应与卷积分析
一、实验目的
加深对离散系统的差分方程、冲激响应与卷积分析方法的理解。
二、实验原理
离散系统
] [n x] [n yDiscrete-timesystme 其输入、输出关系可用以下差分方程描述: MkkNkkk n x p k n y d0 0] [ ] [
输入信号分解为冲激信号, mm n m x n x ] [ ] [ ] [ 。记系统单位冲激响应 ] [ ] [ n h n ,则系统响应为如下的卷积计算式: mm n h m x n h n x n y ] [ ] [ ] [ ] [ ] [
当 N k d k ,... 2 , 1 , 0 时,h[n]就是有限长度的(n:[0,M]),称系统为 FIR 系统;反之,称系统为 IIR 系统。
1、 filter 可用来求一个离散系统的输出。
调用格式:
y=filter(b,a,x); 2、 impz 可用来求一个离散系统的 h(n)。
调用格式:
h=impz(b,a,N);
[h,t]=impz(b,a,N); abnanbz n a z a z az n b z b z b bz Az Bz H ) 1 ( ... ) 3 ( ) 2 ( 1) 1 ( ... ) 3 ( ) 2 ( ) 1 () () () (2 12 1 )] 1 ( , ), 3 ( ), 2 ( ), 1 ( [)] 1 ( , ), 3 ( ), 2 ( ), 1 ( [ abn a a a a an b b b b b 三、实验内容
编制程序求解下列两个系统的输出、单位冲激响应,并绘出其图形。
] 1 [ ] [ ] 2 [ 125 . 0 ] 1 [ 75 . 0 ] [ n x n x n y n y n y; ]} 4 [ ] 3 [ ] 2 [ ] 1 [ { 25 . 0 ] [ n x n x n x n x n y 四、实验要求
给出理论计算结果与程序计算结果并讨论。
1)] 1 [ ] [ ] 2 [ 125 . 0 ] 1 [ 75 . 0 ] [ n x n x n y n y n y
①理论结果: 解: :(1)对于] 1 [ ] [ ] 2 [ 125 . 0 ] 1 [ 75 . 0 ] [ n x n x n y n y n y 计算可得系统函数为:H(z)=2 11125 . 0 75 . 0 1z - 1 z z
=1z 25 . 0 15 -+1z 5 . 0 16 由此可得 h(n)=[-5(-0、25)n+6(-0、5)n]u(n) h(0)=1,
h(1)=-1、75,
h(2)=1、1875,
h(3)=-0、828125 ②程序计算结果: MATLAB 程序如下:
m = [-30:30]; b =[1,-1]; a = [1,0、75,0、125];
x9 = [(m-0)==1];
h =filter(b,a,x9); n = (-10:50);subplot(1,2,1); stem(n,h);axis([-10,50,-1,1、5]); title("Impluse Response"); xlabel("n"),ylabel("h(n)"); subplot(1,2,2);y1 = impz(b,a,n);stem(n,y1); 程序 运行结果 :
图 3 输出函数、单位冲激响应 2) ]} 4 [ ] 3 [ ] 2 [ ] 1 [ { 25 . 0 ] [ n x n x n x n x n y ① 理论结果 : (2)对于]} 4 [ ] 3 [ ] 2 [ ] 1 [ { 25 . 0 ] [ n x n x n x n x n y 由差分方程计算系统函数可得 H(z)=0、25(-4 -3 -2 -1Z Z Z Z ) 则单位冲击响应为:h(n)=0、25(δ(n-1)+δ(n-2)+δ(n-3)+δ(n-4)),该波形与仿真结果一致。
②程序计算结果: MATLAB 程序如下:
m = [-30:30]; b =[0、25,0、25,0、25,0、25]; a =1; x9 = [(m-0)==1];
h =filter(b,a,x9); n = (-10:50); subplot(1,2,1); stem(n,h);axis([-10,50,-1,1、5]); title("Impluse Response"); xlabel("n"),ylabel("h(n)"); subplot(1,2,2);y1 = impz(b,a,n);stem(n,y1);
③程序 运行结果:
图 4 输出函数、单位冲激响应 实验三、离散系统的零、极点分布及频率 响应分析
一、实验目的
加深对离散系统的频率响应分析与零、极点分布的概念理解。
二、实验原理
离散系统的时域方程为
z freqz 可用来求一个离散系统的频率响应。
调用格式: :
[H,w] = f reqz(B B, ,A A ,N,’whole’);
B 与 A 分别为离散系统的(系统函数分子、分母多项式的系数向量); N 表示选取单位圆的上半圆等间距的 N 个点作为频响输出;返回量 H 则包含了离散系统频响在 0—pi 范围内 N 个频率等分点的值,向量 w 则包含范围内 N 个频率等分点。调用中若 N 默认,默认值为 512。
三、实验内容
已知系统差分方程如下: y(n)-1、6y(n-1)+1、28y(n-2) =0、5x(n)+0、1x(n-1) 四、实验要求
(1)编程得到系统频响的幅度响应与相位响应,并画图。
(2)编程得到系统的零极点分布图,分析系统的因果性与稳定性。
(3)给出理论计算结果与程序计算结果并讨论。
提示:幅度响应 mag=abs(H),plot(w,mag) 相位响应 ph=angle(H),plot(w,ph) 五. . 实验结果 分析处理
y(n)-1、6y(n-1)+1、28y(n-2) =0、5x(n)+0、1x(n-1) 1)编程得到系统频响的幅度响应与相位响应,并画图。
①B MATLAB 程序 如下: n=-30:30;b = [0、5,0、1], a =[1,-1、6,1、28]; [H,w] = freqz(b,a,512,"whole"); mag = abs(H),subplot(1,2,1); plot(w,mag);ph = angle(H),subplot(1,2,2);plot(w,ph);
②B MATLAB 运行结果:
图 1
系统频响的幅度响应与相位响应 (2)编程得到系统的零极点分布图,分析系统的因果性与稳定性。
①B MATLAB 程序 如下:
n=-30:30;b = [0、5,0、1],
a =[1,-1、6,1、28];[z,p,K] = tf2zp(b,a);zplane(z,p);
②B MATLAB 运行结果:
图 2
系统的零极点分布图
③理论计算 分析: 解: :由理论计算可得 H(z)=2 - 1128 . 1 z 6 . 1 1z 1 . 0 5 . 0z , 解得:极点为:z1=0、8+0、8j,
z2=0、8-0、8j。零点为:z=-0、2 其中|z|=|z1|=|z2|=0、8 2 =1、1312,当|z|>=0、8 2 时,收敛域包括无穷大,系统就是因果系统;单位圆不在敛域内面,所以系统不稳定。当|z|<0、8 2 时,收敛域不包括无穷大,系统不就是因果系统;单位圆在收敛域内面,所以系统稳定。
实验四、T DFT 算法的应用
一、实验目的
利用 DFT 对信号(如由多个正弦信号组成的信号)进行频谱分析,并研究采样长度、截断(即加窗)与补零对分析频率的影响。
二、实验原理
1. 截断:通常情况下,信号都就是无限长的。而在运用计算机进行模拟时,这就是无法操作的。所以实际情况下,要把观测的信号限制在一定长的时间之内。为了从无限长的信号中得到有限长的数据,在时域乘一个窗函数,将信号截短,叫做加窗。
2. 补零:为了增加频域抽样点数 N,在不改变时域数据的情况下,在时域数据末端加一些零值点,叫做补零。
3. 频率分辨率:指对两个最近的频谱峰值能够分辨的能力。
利用 DFT 进行频谱分析,并研究不同数据长度,补零,加窗等对频率分辨率的影响。
利用 DFT 计算频谱的目的在于,针对计算机只能计算有限个离散的点的取值这一特点,实现计算机对连续时间信号的频谱的模拟。所以我们比较关心的就是模拟频谱与原信号频谱的拟合程度,我们希望拟合程度越高越好。这就需要增加频率分辨率,因为频率分辨率越高,根据公式,说明相同采样频率下,采样的长度就越长,也就就是频谱采样的点数就越多,我们可以瞧到的模拟频谱图像就越清晰,这样与原信号的拟合程度就越好。
MATLAB 提供了 2 个内部函数用于计算 DFT,它们分别就是:
fft(x,N)
fft(x,N) 计算 N 点的 DFT。
三、实验内容
已知序列 x=sin(5*2*pi*n/fs)+cos(3*2*pi*n/fs),fs=30 (1)不同采样长度的影响:分别计算序列的 32 点与 128 点 DFT,绘出幅度谱图形。
F32=fft(x,32); F32mag=abs(F32); (2)补零对频率分辨率的影响:设定序列长度为 32,采样长度 64,对连续时间信号进行采样,首先对采样的信号时域图像后补若干个零,然后计算出采样长度情况下的频谱,并给出 32 点与 64 点的序列图与幅度谱图。
m1=64;m=0:m1-1 x2=[x1,zeros(1,m1-n1)] (3)加窗对频率分辨率的影响:一部分就是设定采样长度 N1,对连续时间信号进行采样,然后利用 DFT 计算出采样长度 N1 情况下的频谱,并给出时域与频域图像。一部分采用加窗的方法对时域图像进行采样,然后利用 DFT 计算出采样长度 N1 情况下的频谱,并给出时域与频域图像。这里要注意的就是我们采用加窗的方法时,设定信号的长度为 400,然后对长度为 400 的信号进行截短加窗,这时只有加窗处我们就是对信号采样的,即得到的信号就是有效的,其余位置信号都为0,所以这时我们仍然认为采样长度为 128、 加窗: n=0:399;m=0:399;m2=128; fs=30; x= sin(5*2*pi*n/fs)+cos(3*2*pi*n/fs); x1=x、*(heaviside(m)-heaviside(m-m2)); F1=fft(x1); mag=abs(F1); stem(mag);
四. . 实验结果 分析处理
已知序列 x=sin(5*2*pi*n/fs)+cos(3*2*pi*n/fs),fs=30 (1)不同采样长度的影响:分别计算序列的 32 点与 128 点 DFT,绘出幅度谱图形。
①B MATLAB 程序 如下: fs = 30; n = 0:31; x1 = sin(5*2*pi*n/fs) + cos(3*2*pi*n/fs); F32 = fft(x1,32);F32mag = abs(F32); subplot(1,2,1);stem(n,F32mag);
m = 0:127;x2 = sin(5*2*pi*m/fs) + cos(3*2*pi*m/fs); F128 = fft(x2,128); F128mag = abs(F128); subplot(1,2,2); stem(m,F128mag); ②B MATLAB 运行结果:
图 1
32 点与 128 点 DFT 幅度谱图形
(2)补零对频率分辨率的影响:设定序列长度为 32,采样长度 64,对连续时间信号进行采样,首先对采样的信号时域图像后补若干个零,然后计算出采样长度情况下的频谱,并给出 32 点与 64 点的序列图与幅度谱图。
①B MATLAB 程序 如下: fs = 30;n = 0:31; x1 = sin(5*2*pi*n/fs) + cos(3*2*pi*n/fs); F32 = fft(x1,32); F32mag = abs(F32); subplot(1,2,1);stem(n,x1);
subplot(1,2,2);stem(n,F32mag); m1= 64;n1 = 32; x2=[x1,zeros(1,m1-n1)]; F64 = fft(x2,64); F64mag = abs(F64); subplot(1,2,1);stem(n2,x2); subplot(1,2,2);stem(n2,F64mag); ②B MATLAB 运行结果:
图 2 补零后 32 点与 64 点的序列图与幅度谱图 (3)加窗对频率分辨率的影响:一部分就是设定采样长度 N1,对连续时间信号进行采样,然后利用 DFT 计算出采样长度 N1 情况下的频谱,并给出时域与频域图像。一部分采用加窗的方法对时域图像进行采样,然后利用 DFT 计算出采样长度N1 情况下的频谱,并给出时域与频域图像。
①B MATLAB 程序 如下:
fs = 30;n = 0:399;m = 0:399;m2=128;x3 = sin(5*2*pi*n/fs) + cos(3*2*pi*n/fs); x5=[ones(1,128),zeros(1,272)]; x4 = x3、*x5; F1 = fft(x4); mag = abs(F1);stem(mag); ②B MATLAB 运行结果:
图 3
400 点 DFT 序列图 分析:由结果图可知,加穿函数后信号的频谱图近似等于信号原谱图。
实验五R :IIR 滤波器的设计
一、实验目的
1.理解滤波器参数的意义; 2.掌握脉冲响应不变法与双线性变换法设计 IIR 数字滤波器的方法; 3.掌握利用 Matlab 设计其它各型 IIR 数字滤波器的方法; 4.掌握分析滤波器就是否达到性能指标的方法。
二、实验原理
1、利用脉冲响应不变法,直接根据归一化的巴特沃斯低通模拟滤波器系统函数H ( p )得到 IIR 数字低通滤波器方法就是
R IIR 数字滤波器的设计步骤: :
(1)按照一定规则把给定的滤波器技术指标转换为模拟低通滤波器的技术指标;
模拟通带、阻带截止频率
(2)根据模拟滤波器技术指标设计为响应的模拟低通滤波器;
[n,Wn]=buttord(wp,ws,Ap,As,’s’) [z,p,k]=buttap(n);设计模拟滤波器 [b,a]=zp2tf(z,p,k);
将系统函数的零极点转化为系统函数一般形式的系数 [bn,an]=lp2lp(b,a,Wn); (3)根据脉冲响应不变法与双线性不变法把模拟滤波器转换为数字滤波器; [bz,az]=bilinear(bn,an,1/Ts) [bz,az]=
impinvar (bn,an,1/Ts)
(4)如果要设计的滤波器就是高通、带通或带阻滤波器,则首先把它们的技术指标转化为模拟低通滤波器的技术指标,设计为数字低通滤波器,最后通过频率转换的方法来得到所要的滤波器。
3、直接设计数字滤波器方法: [n,Wn]=buttord(wp/pi,ws/pi,Ap,As) [b,a]=butter(n,Wn,’high’) 设计其它各型 IIR 数字滤波器的理论方法在这里不再给出,读者可参瞧有关内容。在 Matlab 中,设计滤波器的过程很简单,只要加上一些控制字符即可。控制字符省略或为“low”表示设计低通滤波器,控制字符为“high”表示设计高通滤波器,控制字符为“band”表示设计带通滤波器,控制字符为“stop”表示设计带阻滤波器。
三、实验内容
freqz(b,a); figure; [y,t]=impz(bz,az,101); stem(t,y) 四、实验 结果
(1)脉冲响应不变法设计:
①
B MATLAB 程序 : Ts = 0、001; [n1,Wn1] = buttord(0、2*pi/Ts,0、6*pi/Ts,2,15,"s"), [z,p,k] = buttap(n1),
[b,a] = zp2tf(z,p,k); [b1,a1] = lp2lp(b,a,Wn1); [bz1,az1] = impinvar(b1,a1,1000); freqz(bz1,az1);
figure; [y1,t1] = impz(bz1,az1,101); stem(t1,y1); ②
B MATLAB 仿真结果 :
图 1
As=15dB 频率响应截图
图 2
冲击响应图 ③
理论 计算 分析 :
ωp = 0、2pi 时有
Ap =< 2dB ;ωs = 0、6pi 时有
As > 15dB,符合设计要求。
(2)双线性变换法设计
①
B MATLAB 程序 : Ts = 0、001;
[n1,Wn1] = buttord(0、2*pi/Ts,0、6*pi/Ts,2,15,"s"), [z,p,k] = buttap(n1), [b,a] = zp2tf(z,p,k); [b1,a1] = lp2lp(b,a,Wn1); [bz1,az1] = bilinear(b1,a1,1000); freqz(bz1,az1); figure; [y1,t1] = impz(bz1,az1,101); stem(t1,y1); ②
B MATLAB 仿真结果 :
图 3
As=18dB 频率响应截图
图 4
冲击响应图 ③
理论分析 :
ωp = 0、2pi 时有
Ap =< 2dB ;ωs = 0、6pi 时有
As > 15dB,符合设计要求。
(3)直接设计法设计高通滤波器 ①
B MATLAB 程序 : Ts = 0、001; [n1,Wn1] = buttord(0、8,0、44,3,20), [b,a] = butter(n1,Wn1,"high"); freqz(b,a); figure; [y1,t1] = impz(b,a,101); stem(t1,y1); ②
B MATLAB 仿真结果: :
图 6
As=20dB 频率响应截图
图 7
冲击响应截图
③
理 论分析 :
ωp = 0、8pi 时有
Ap =< 3dB ;ωs = 0、44pi 时有
As > 20dB,符合设计要求。
实验六、R FIR 滤波器的设计
一、实验目的
1．掌握 FIR 数字滤波器的设计方法; 2．熟悉利用 MATLAB 软件进行 FIR 数字滤波器设计,以及对所设计的滤波器进行分析; 二、 实验原理
窗函数法的设计步骤
窗函数设计 R FIR 数字滤波器就是从时域出发的, , 把理想滤波器的单位取样响应用于合适的窗函数截短成为有限长度的 h(n), 使得 h(n) 逼近理想的 h h d d (n) 。以实现所设计的滤波器的频率响应 H(ejw) ) 逼近与理想滤波器的频率响应 Hd(ejw) 。
由过渡带宽度及阻带最小衰减要求, , 可选定窗形状, , 并估计窗口长度 N N 。
b Matlab 提供了 1 fir1 函数, , 以实现线性相位 R FIR 滤波器。
调用格式如下: :
b=fir1(N,wc, ‘ ftype ’ ,window(M))
参数: :
N: 阶数
wc: 归一化的数字频率 ,0 ≤ wc ≤ 1
wc=(wp+(ws- - wp)/2)/pi
ftype: 滤波器类型, , 如高通、带阻等。
window: 应用的窗函数类型
实验内容
基于窗函数设计法,利用 MATLAB 软件设计满足设计要求的 FIR 数字低通滤波器。
(1)要求利用窗函数设计法设计 FIR 数字低通滤波器,滤波器参数要求为:
。要求设计数字低通滤波器,同时要求给出滤波器的辐频特性曲线图以及幅度响应图与对数幅度响应图; 参考函数: freqz(b,a); [h,w]=freqz(b,a); plot(w,abs(h)) plot(w,20*log10(abs(h))) 三、实验结果 （1）
汉宁窗
①B MATLAB 程序: wp = 0、2*pi;ws = 0、3*pi;a = 1; wc = (wp+(ws-wp)/2)/pi; b = fir1(62,wc,"low",hanning(63)); freqz(b,a); figure; [h,w] = freqz(b,a); plot(w,abs(h)); figure; plot(w,20*log10(abs(h)));
②B MATLAB 运行 结果:
图 1
辐频特性曲线图
图 2
幅度响应图
图 3 对数幅度响应图 ③理论分析计算: 当 wp=0、2*pi 时,Ap=<0、25dB,当 ws=0、3*pi,As=>50dB。阻带衰减符合要求。
（2）
海明窗
①B MATLAB 程序: wp = 0、2*pi;ws = 0、3*pi;a = 1; wc = (wp+(ws-wp)/2)/pi; b = fir1(66,wc,"low",hamming(67)); freqz(b,a); figure; [h,w] = freqz(b,a); plot(w,abs(h)); figure; plot(w,20*log10(abs(h)));
②B MATLAB 运行 结果:
图 4
辐频特性曲线图
图 5
幅度响应图
图 6 对数幅度响应图 ③理论分析计算: :
当 wp=0、2*pi 时,Ap=<0、25dB,当 ws=0、3*pi,As<50dB。阻带衰减稍微不符合要求。
（3）
布莱克曼窗
①B MATLAB 程序: wp = 0、2*pi;ws = 0、3*pi;a = 1; wc = (wp+(ws-wp)/2)/pi; b = fir1(110,wc,"low",Blackman(111)); freqz(b,a); figure; [h,w] = freqz(b,a); plot(w,abs(h)); figure; plot(w,20*log10(abs(h)));
②B MATLAB 运行 结果:
图 7
辐频特性曲线图
图 8
幅度响应图
图 9 对数幅度响应图 ③ 理论分析计算: 当 wp=0、2*pi 时,Ap=<0、25dB,当 ws=0、3*pi,As=>50dB。阻带衰减符合要求,甚至远大于要求衰减。
四、实验分析
窗函数的选择原则就是保证阻带衰减的情况下选择主瓣窄的窗函数。由于滤波器的参数要求就是,海明窗的阻带最小衰减为-53db,满足该滤波器阻带衰减要求,海宁窗为-44db,不满足要求,布拉克曼窗为-74db 也满足要求,理论上说海明窗最适合设计该滤波器。由程序结果可知,海宁窗与海明窗的主瓣宽度一致,而布莱克曼窗主瓣宽度比她们大,但就是阻带衰减最大。
上一篇：热敏电阻实验报告