数字信号处理实验.docx
- 文档编号:26142148
- 上传时间:2023-06-17
- 格式:DOCX
- 页数:26
- 大小:758.80KB
数字信号处理实验.docx
《数字信号处理实验.docx》由会员分享,可在线阅读,更多相关《数字信号处理实验.docx(26页珍藏版)》请在冰豆网上搜索。
数字信号处理实验
太原理工大学
数字信号处理课程实验报告
专业班级
学号2013000000
姓名XXX
指导教师XXX
实验一:
系统响应及系统稳定性
1.实验目的
(1)掌握 求系统响应的方法。
(2)掌握时域离散系统的时域特性。
(3)分析、观察及检验系统的稳定性。
2.实验原理与方法
在时域中,描写系统特性的方法是差分方程和单位脉冲响应,在频域可以用系统函数描述系统特性。
已知输入信号可以由差分方程、单位脉冲响应或系统函数求出系统对于该输入信号的响应,本实验仅在时域求解。
在计算机上适合用递推法求差分方程的解,最简单的方法是采用MATLAB语言的工具箱函数filter函数。
也可以用MATLAB语言的工具箱函数conv函数计算输入信号和系统的单位脉冲响应的线性卷积,求出系统的响应。
系统的时域特性指的是系统的线性时不变性质、因果性和稳定性。
重点分析实验系统的稳定性,包括观察系统的暂态响应和稳定响应。
系统的稳定性是指对任意有界的输入信号,系统都能得到有界的系统响应。
或者系统的单位脉冲响应满足绝对可和的条件。
系统的稳定性由其差分方程的系数决定。
实际中检查系统是否稳定,不可能检查系统对所有有界的输入信号,输出是否都是有界输出,或者检查系统的单位脉冲响应满足绝对可和的条件。
可行的方法是在系统的输入端加入单位阶跃序列,如果系统的输出趋近一个常数(包括零),就可以断定系统是稳定的[19]。
系统的稳态输出是指当时,系统的输出。
如果系统稳定,信号加入系统后,系统输出的开始一段称为暂态效应,随n的加大,幅度趋于稳定,达到稳态输出。
注意在以下实验中均假设系统的初始状态为零。
3.实验内容及步骤
(1)编制程序,包括产生输入信号、单位脉冲响应序列的子程序,用filter函数或conv函数求解系统输出响应的主程序。
程序中要有绘制信号波形的功能。
(2)给定一个低通滤波器的差分方程为y(n)=0.05*x(n)+0.05*x(n-1)+0.9*y(n-1)
输入信号x1(n)=R8(n),x2(n)=u(n)
a)分别求出x1(n)=R8(n)和x2(n)=u(n)的系统响应,并画出其波形。
b)求出系统的单位冲响应,画出其波形。
(3)给定系统的单位脉冲响应为
h1(n)=R10(n)
h2(n)=δ(n)+2.5*δ(n-1)+2.5*δ(n-2)+δ(n-3)
用线性卷积法求x1(n)=R8(n)分别对系统h1(n)和h2(n)的输出响应,并画出波形。
(4)给定一谐振器的差分方程为
y(n)=1.8237*y(n-1)-0.9801*y(n-2)+b0*x(n)-b0*x(n-2)
令b0=1/100.49 ,谐振器的谐振频率为0.4rad。
a)用实验方法检查系统是否稳定。
输入信号为u(n)时,画出系统输出波形。
b)给定输入信号为x(n)=sin(0.014*n)+sin(0.4*n),求出系统的输出响应,并画出其波形。
4.实验程序清单以及波形图
(2)y(n)=0.05x(n)+0.05x(n-1)+0.9y(n-1)输入信号为x1(n)=R8(n),x2(n)=u(n)
B=[0.05,0.05];
A=[1,-0.9];
x1n=[11111111zeros(1,100)];
x2n=[ones(1,100)];
hn=impz(B,A,58);
subplot(2,2,1);
stem(hn,'r','.');
xlabel('n');ylabel('h(n)');
title('(a)系统单位冲击响应');
y1n=filter(B,A,x1n);
subplot(2,2,2);
stem(y1n,'g','.');
xlabel('n');ylabel('y1n');
title('(b)系统对R8(n)的响应y(n)');
y2n=filter(B,A,x2n);
subplot(2,2,3);
stem(y2n,'y','.');
xlabel('n');ylabel('y2n');
title('(c)系统对a(n)的响应y2(n)');
(3)h1(n)=R10(n),h2n=δ(n)+2.5δ(n-1)+2.5δ(n-2)+δ(n-3),x1(n)=R8(n)
(4)y(n)=1.8237y(n-1)-0.9801y(n-2)+b0x(n)-b0x(n-2),b0=1/100.49
谐振频率0.4rad,输入信号x(n)=sin(0.014n)+sin(0.4n)
5.思考题
(1)如果输入信号为无限长序列,系统的单位脉冲响应是有限长序列,可否用线性卷积法求系统的响应?
如何求?
答:
对输入信号序列分段;求单位脉冲响应h(n)与各段的卷积;将各段卷积结果相加。
具体实现方法有重叠相加法和重叠保留法。
(2)如果信号经过低通滤波器,把信号的高频分量滤掉,时域信号会有何变化?
用前面的第一个实验结果进行分析说明。
答:
如果信号经过低通滤波器,把信号的高频分量滤掉,时域信号将被平滑。
由实验
(1)的实验图可见,经过系统低通滤波器滤波使输入信号δ(n),x1(n)=R8(n)和x2(n)=u(n)的阶越变化变得缓慢上升与下降。
实验三:
用FFT对信号作频谱分析
1.实验目的
学习用FFT对连续信号和时域离散信号进行谱分析的方法,了解可能出现的分析 误差及其原因,以便正确应用FFT。
2.实验原理
用FFT对信号作频谱分析是学习数字信号处理的重要内容。
经常需要进行谱分析的信号是模拟信号和时域离散信号。
对信号进行谱分析的重要问题是频谱分辨率D和分析误差。
频谱分辨率直接和FFT的变换区间N有关,因为FFT能够实现的频率分辨率是2pi/N,因此要求2pi/N≦D。
可以根据此式选择FFT的变换区间N。
误差主要来自于用FFT作频谱分析时,得到的是离散谱,而信号(周期信号除外)是连续谱,只有当N较大时离散谱的包络才能逼近于连续谱,因此N要适当选择大一些。
周期信号的频谱是离散谱,只有用整数倍周期的长度作FFT,得到的离散谱才能代表周期信号的频谱。
如果不知道信号周期,可以尽量选择信号的观察时间长一些。
对模拟信号进行谱分析时,首先要按照采样定理将其变成时域离散信号。
如果是模拟周期信号,也应该选取整数倍周期的长度,经过采样后形成周期序列,按照周期序列的谱分析进行。
3.实验步骤及内容
(1)对以下序列进行谱分析。
X1(n)=R4(n)
X2(n)={n+10≦n≦3
8-n4≦n≦7
0其它n}
X3(n)={4-n0≦n≦3
n-34≦n≦7
0其它n}
选择FFT的变换区间N为8和16两种情况进行频谱分析。
分别打印其幅频特性曲线。
并进行对比、分析和讨论。
(2)对以下周期序列进行谱分析:
X4(n)=cos(pi/4)n
X5(n)=cos(pi/4)n+cos(pi/8)n
选择FFT的变换区间N为8和16两种情况分别对以上序列进行频谱分析。
分别打印其幅频特性曲线。
并进行对比、分析和讨论。
(3)对模拟周期信号进行谱分析:
X8(t)=cos(8*pi*t)+cos(16*pi*t)+cos(20*pi*t)
选择采样频率Fs=64Hz,对变换区间N=16,32,64三种情况进行谱分析。
分别打印其幅频特性,并进行分析和讨论。
4.实验程序清单以及波形图
(1)x1n=[1,1,1,1];
M=8;xa=1:
(M/2);
xb=(M/2):
-1:
1;
x2n=[xa,xb];
x3n=[xb,xa];
x1k8=fft(x1n,8);x1k16=fft(x1n,16);
x2k8=fft(x2n,8);x2k16=fft(x2n,16);
x3k8=fft(x3n,8);x3k16=fft(x3n,16);
n=0:
length(x1k8)-1;k1=2*n/length(x1k8);
subplot(3,2,1);stem(k1,abs(x1k8),'.');
title('(1a)8点DFT[x_1(n)]');
xlabel('w/pi');ylabel('幅度');
axis([0,2,0,1.2*max(abs(x1k8))]);
n=0:
length(x1k16)-1;k2=2*n/length(x1k16);
subplot(3,2,2);stem(k2,abs(x1k16),'.');
title('(1b)16点DFT[x_1(n)]');
xlabel('ω/pi');ylabel('幅度');
axis([0,2,0,1.2*max(abs(x1k16))])
n=0:
length(x2k8)-1;k3=2*n/length(x2k8);
subplot(3,2,3);stem(k3,abs(x2k8),'.');
title('(2a)8点DFT[x_2(n)]');
xlabel('ω/pi');ylabel('幅度');
axis([0,2,0,1.2*max(abs(x2k8))])
n=0:
length(x2k16)-1;k4=2*n/length(x2k16);
subplot(3,2,4);stem(k4,abs(x2k16),'.');
title('(2b)16点DFT[x_2(n)]');
xlabel('ω/pi');ylabel('幅度');
axis([0,2,0,1.2*max(abs(x2k16))])
n=0:
length(x3k8)-1;k5=2*n/length(x3k8);
subplot(3,2,5);stem(k5,abs(x3k8),'.');
title('(3a)8点DFT[x_3(n)]')
;xlabel('ω/pi');ylabel('幅度');
axis([0,2,0,1.2*max(abs(x3k8))])
n=0:
length(x3k16)-1;k6=2*n/length(x3k16);
subplot(3,2,6);stem(k6,abs(x3k16),'.');
title('(3b)16点DFT[x_3(n)]');
xlabel('ω/pi');ylabel('幅度');
axis([0,2,0,1.2*max(abs(x3k16))])
(2)N=8;n=0:
N-1;
x4n=cos(pi*n/4);
x5n=cos(pi*n/4)+cos(pi*n/8);
X4k8=fft(x4n);
X5k8=fft(x5n);
N=16;n=0:
N-1;
x4n=cos(pi*n/4);
x5n=cos(pi*n/4)+cos(pi*n/8);
X4k16=fft(x4n);
X5k16=fft(x5n);
b=0:
length(X4k8)-1;l1=2*b/length(X4k8);
subplot(2,2,1);stem(l1,X4k8,'.');
title('(4a)8点DFT[x_4(n)]');
xlabel('ω/pi');ylabel('幅度');
axis([0,2,0,1.2*max(abs(X4k8))])
b=0:
length(X4k16)-1;l2=2*b/length(X4k16);
subplot(2,2,2);stem(l2,X4k16,'.');
title('(4b)16点DFT[x_4(n)]');
xlabel('ω/pi');ylabel('幅度');
axis([0,2,0,1.2*max(abs(X4k16))])
b=0:
length(X5k8)-1;l3=2*b/length(X5k8);
subplot(2,2,3);stem(l3,X5k8,'.');
title('(5a)8点DFT[x_5(n)]');
xlabel('ω/pi');ylabel('幅度');
axis([0,2,0,1.2*max(abs(X5k8))])
b=0:
length(X5k16)-1;l4=2*b/length(X5k16);
subplot(2,2,4);stem(l4,X5k16,'.');
title('(5b)1点DFT[x_5(n)]');
xlabel('ω/pi');ylabel('幅度');
axis([0,2,0,1.2*max(abs(X5k16))])
(3)figure(4)
Fs=64;T=1/Fs;
N=16;n=0:
N-1;
x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);
X6k16=fft(x6nT);
X6k16=fftshift(X6k16);
Tp=N*T;F=1/Tp;
k=-N/2:
N/2-1;fk=k*F;
subplot(3,1,1);stem(fk,abs(X6k16),'.');boxon
title('(6a)16点|DFT[x_6(nT)]|');xlabel('f(Hz)');ylabel('幅度');
axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k16))])
N=32;n=0:
N-1;
x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);
X6k32=fft(x6nT);
X6k32=fftshift(X6k32);
Tp=N*T;F=1/Tp;
k=-N/2:
N/2-1;fk=k*F;
subplot(3,1,2);stem(fk,abs(X6k32),'.');boxon
title('(6b)32点|DFT[x_6(nT)]|');xlabel('f(Hz)');ylabel('幅度');
axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k32))])
N=64;n=0:
N-1;
x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T);
X6k64=fft(x6nT);
X6k64=fftshift(X6k64);
Tp=N*T;F=1/Tp;
k=-N/2:
N/2-1;fk=k*F;
subplot(3,1,3);stem(fk,abs(X6k64),'.');boxon
title('(6a)64点|DFT[x_6(nT)]|');xlabel('f(Hz)');ylabel('幅度');
axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k64))])
5.思考题
(1)对于周期序列,如果周期不知道,如何用FFT进行谱分析?
答:
周期信号的周期预先不知道时,可先截取M点进行DFT,再将截取长度扩大1倍截取,比较结果,如果二者的差别满足分析误差要求,则可以近似表示该信号的频谱,如果不满足误差要求就继续将截取长度加倍,重复比较,直到结果满足要求。
(2)如何选择FFT的变换区间?
(包括非周期信号和周期信号)
答:
对于非周期信号:
有频谱分辨率F,而频谱分辨率直接和FFT的变换区间有关,因为FFT能够实现的频率分辨率是2π/N...因此有最小的N>2π/F。
就可以根据此式选择FFT的变换区间。
对于周期信号,周期信号的频谱是离散谱,只有用整数倍周期的长度作FFT,得到的离散谱才能代表周期信号的频谱。
(3)当N=8时,和的幅频特性会相同吗?
为什么?
N=16呢?
答:
不相同,由DFT的定义可知,x(n)的作用是对相同的K值的幅度值,所以虽然不影响频谱分布,但是影响频率幅值.由于xn2和xn3的数值分布范围不同,所以导致在相同K值时幅值不同,因此xn2和xn3的幅频特性不相同.N=16时也不相同。
实验四:
IIR数字滤波器设计及软件实现
1.实验目的
(1)熟悉用双线性变换法设计IIR数字滤波器的原理与方法;
(2)学会调用MATLAB信号处理工具箱中滤波器设计函数(或滤波器设计分析工具fdatool)设计各种IIR数字滤波器,学会根据滤波需求确定滤波器指标参数。
(3)掌握IIR数字滤波器的MATLAB实现方法。
(4)通过观察滤波器输入输出信号的时域波形及其频谱,建立数字滤波的概念。
2.实验原理
设计IIR数字滤波器一般采用间接法(脉冲响应不变法和双线性变换法),应用最广泛的是双线性变换法。
基本设计过程是:
①先将给定的数字滤波器的指标转换成过渡模拟滤波器的指标;②设计过渡模拟滤波器;③将过渡模拟滤波器系统函数转换成数字滤波器的系统函数。
MATLAB信号处理工具箱中的各种IIR数字滤波器设计函数都是采用双线性变换法。
第六章介绍的滤波器设计函数butter、cheby1、cheby2和ellip可以分别被调用来直接设计巴特沃斯、切比雪夫1、切比雪夫2和椭圆模拟和数字滤波器。
本实验要求读者调用如上函数直接设计IIR数字滤波器。
本实验的数字滤波器的MATLAB实现是指调用MATLAB信号处理工具箱函数filter对给定的输入信号x(n)进行滤波,得到滤波后的输出信号y(n)。
3.实验内容及步骤
(1)调用信号产生函数mstg产生由三路抑制载波调幅信号相加构成的复合信号st,该函数还会自动绘图显示st的时域波形和幅频特性曲线,如图10.4.1所示。
由图可见,三路信号时域混叠无法在时域分离。
但频域是分离的,所以可以通过滤波的方法在频域分离,这就是本实验的目的。
图10.4.1 三路调幅信号st的时域波形和幅频特性曲线
(2)要求将st中三路调幅信号分离,通过观察st的幅频特性曲线,分别确定可以分离st中三路抑制载波单频调幅信号的三个滤波器(低通滤波器、带通滤波器、高通滤波器)的通带截止频率和阻带截止频率。
要求滤波器的通带最大衰减为0.1dB,阻带最小衰减为60dB。
提示:
抑制载波单频调幅信号的数学表示式为S(t)=cos(2*pi*f0*t)*cos(2*pi*fc*t)=1/2*[cos(2*pi(fc-f0)*t)+cos(2*pi(fc+f0)*t)
其中,cos(2*pi*fc*t)称为载波,fc为载波频率,cos(2*pi*f0*t)称为单频调制信号,f0为调制正弦波信号频率,且满足fc>f0。
由上式可见,所谓抑制载波单频调幅信号,就是2个正弦信号相乘,它有2个频率成分:
和频fc+f0和差频fc-f0,这2个频率成分关于载波频率fc对称。
所以,1路抑制载波单频调幅信号的频谱图是关于载波频率fc对称的2根谱线。
容易看出,图10.4.1中三路调幅信号的载波频率分别为250Hz、500Hz、1000Hz。
有关调幅(AM)和抑制载波调幅(SCAM)的一般原理与理念,请参考通信原理教材。
(3)编程序调用MATLAB滤波器设计函数ellipord和ellip分别设计这三个椭圆滤波器,并绘图显示其幅频响应特性曲线。
(4)调用滤波器实现函数filter,用三个滤波器分别对信号产生函数mstg产生的信号st进行滤波,分离出st中的三路不同载波频率的调幅信号y1(n)、y2(n)和y3(n),并绘图显示y1(n)、y2(n)和y3(n)的时域波形,观察分离效果。
4.信号产生函数mstg清单
functionst=mstg
%产生信号序列向量st,并显示st的时域波形和频谱
%st=mstg返回三路调幅信号相加形成的混合信号,长度N=800
N=800 %N为信号st的长度。
Fs=10000;T=1/Fs;Tp=N*T;%采样频率Fs=10kHz,Tp为采样时间
t=0:
T:
(N-1)*T;k=0:
N-1;f=k/Tp;
fc1=Fs/10; %第1路调幅信号的载波频率fc1=1000Hz,
fm1=fc1/10; %第1路调幅信号的调制信号频率fm1=100Hz
fc2=Fs/20; %第2路调幅信号的载波频率fc2=500Hz
fm2=fc2/10; %第2路调幅信号的调制信号频率fm2=50Hz
fc3=Fs/40; %第3路调幅信号的载波频率fc3=250Hz,
fm3=fc3/10; %第3路调幅信号的调制信号频率fm3=25Hz
xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t);%产生第1路调幅信号
xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t);%产生第2路调幅信号xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t);%产生第3路调幅信号
st=xt1+xt2+xt3; %三路调幅信号相加
fxt=fft(st,N); %计算信号st的频谱
%====以下为绘图部分,绘制st的时域波形和幅频特性曲线====================
subplot(3,1,1)
plot(t,st);grid;xlabel('t/s');ylabel('s(t)');
axis([0,Tp/8,min(st),max(st)]);title('(a)s(t)的波形')
subplot(3,1,2)
stem(f,abs(fxt)/max(abs(fxt)),'.');grid;title('(b)s(t)的频谱')
axis([0,Fs/5,0,1.2]);
xlabel('f/Hz');ylabel('幅度')
5.实验程序框图
调用函数mstg产生st,自动绘图显示st的时域波形图和幅频特性曲线
↓
调用ellipord和ellip分别设计三个椭圆滤波器,并绘图显示其幅频响应特性曲线
↓
调用filter,用三个滤波器分别对信号st进行滤波,分离出三路不同载波频率的调幅信号y1(n),y2(n),和y3(n)
↓
绘图显示y1(n),y2(n),y3(n)的时域波形
↓
结束
6.试验程序清单及波形图
functionst=mstg
N=800
Fs=10000;T=1/Fs;Tp=N*T;
t=0:
T:
(N-1)*T;k=0:
N-1;f=k/Tp;
fc1=Fs/10;
fm1=fc1/10;
fc2=Fs/20;
fm2=fc2/10;
fc3=Fs/40;
fm3=fc3/10;
xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t);
xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t);
xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t);
st=xt1+xt2+xt3;
fxt=fft(st,N);
subplot(3,1,1)
plot(t,st);grid;xlabel('t/s');ylabel('s(t)');
axis([0,Tp/8,min(st),max(st)]);title('(a)s(t)的波形')
su
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数字信号 处理 实验