吉林大学数字信号处理实验报告.docx
《吉林大学数字信号处理实验报告.docx》由会员分享，可在线阅读，更多相关《吉林大学数字信号处理实验报告.docx（39页珍藏版）》请在冰点文库上搜索。
吉林大学数字信号处理实验报告
数字信号处理课程设计实验报告
基础实验
实验一离散时间系统及离散卷积
一、实验目的
〔1〕熟悉MATLAB软件的使用方法。
〔2〕熟悉系统函数的零极点分布、单位脉冲响应和系统频率响应等概念。
〔3〕利用MATLAB绘制系统函数的零极点分布图、系统频率响应和单位脉冲响应。
〔4〕熟悉离散卷积的概念，并利用MATLAB计算离散卷积。
二、实验内容
1、离散时间系统的单位脉冲响应
〔1〕选择一个离散时间系统；
〔2〕用笔进行差分方程的递推计算；
〔3〕编制差分方程的递推计算程序；
〔4〕在电脑上实现递推运算；
〔5〕将程序计算结果与笔算的计算结果进行比较，验证程序运行的正确性；
2.离散系统的幅频、相频的分析方法
〔1〕给定一个系统的差分方程或单位取样响应；
〔2〕用笔计算几个特殊的幅频、相频的值，画出示意曲线图；
〔3〕编制离散系统的幅频、相频的分析程序；
〔4〕在电脑上进行离散系统的幅频、相频特性，并画出曲线；
〔5〕通过比较，验证程序的正确性；
3.离散卷积的计算
〔1〕选择两个有限长序列，用笔计算其线性卷积；
〔2〕编制有限长序列线性卷积程序；
〔3〕利用计算程序对〔1〕选择的有限长序列进行卷积运算；
〔4〕比较结果验证程序的正确性。
三、实验要求
〔1〕编制实验程序，并给编制程序加注释；
〔2〕按照实验内容项要求完成笔算结果；
〔3〕验证编制程序的正确性，记录实验结果。
〔4〕至少要求一个除参考实例以外的实例，在实验报告中，要描述清楚实例中的系统，并对实验结果进行解释说明。
四、实验程序及其结果
1、离散时间系统的脉冲响应
%y（n）-0.9y（n-1）+0.5y（n-2）=x（n）+0.5x（n-1）设x（n）=2^n;0<=n<=15
clearall;
b=[1,0.5];a=[1,-0.9,0.5];
n=0:
15;x=impseq（0,0,15）
h=filter（b,a,x）;%系统冲击响应
subplot（2,1,1）;stem（n,h）;title（'系统冲激响应'）;
xlabel（'n'）;ylabel（'h'）;
n=0:
15;x=2.*n;nx=0:
15;nh=0:
15;
y=conv_m（x,nx,h,nh）;ny=length（y）;
n=0:
ny-1;
subplot（2,1,2）;stem（n,y）;title（'系统对x（n）响应'）;
xlabel（'n'）;ylabel（'y'）;
B=roots（b）;A=roots（a）;
figure;zplane（B,A）;
2、离散系统的幅频、相频的分析方法
%差分方程为y（n）-1.76y（n-1）+1.1829y（n-2）-0.2781y（n-3）
%=0.0181x（n）+0.0543x（n-1）+0.0543x（n-2）+0.0181x（n-3）
b=[0.0181,0.0543,0.0543,0.0181];
a=[1.000,-1.76,1.1829,-0.2781];
m=0:
length（b）-1;
l=0:
length（a）-1;
K=500;
k=0:
1:
K;
w=pi*k/K;
H=（b*exp（-j*m'*w））./（a*exp（-j*l'*w））;%计算频率相应
magH=abs（H）;%magH为幅度
angH=angle（H）;%angH为相位
subplot（2,1,1）;plot（w/pi,magH）;grid;
xlabel（'以\pi为单位的频率'）;ylabel（'幅度'）;
title（'幅度响应'）;
subplot（2,1,2）;plot（w/pi,angH）;grid;
xlabel（'以\pi为单位的频率'）;ylabel（'相位'）;
title（'相位响应'）;
3、离散卷积的计算
%x=[1,4,3,5,3,6,5],-4<=n<=2
%h=[3,2,4,1,5,3],-2<=n<=3
%求两序列的卷积
clearall;
x=[1,4,3,5,3,6,5];nx=-4:
2;
h=[3,2,4,1,5,3];nh=-2:
3;
ny=（nx
（1）+nh
（1））:
（nx（length（x））+nh（length（h）））;
y=conv（x,h）;
n=length（ny）;
subplot（3,1,1）;stem（nx,x）;xlabel（'nx'）;ylabel（'x'）;
subplot（3,1,2）;stem（nh,h）;xlabel（'nh'）;ylabel（'h'）;
subplot（3,1,3）;stem（ny,y）;xlabel（'n'）;ylabel（'x和h的卷积'）;
实验二离散傅里叶变换与快速傅里叶变换
一、实验目的
〔1〕加深理解离散傅立叶变换及快速傅立叶变换概念；
〔2〕学会应用FFT对典型信号进行频谱分析的方法；
〔3〕研究如何利用FFT程序分析确定性时间连续信号；
〔4〕熟悉应用FFT实现两个序列的线性卷积的方法。
二、实验内容
1、用离散傅立叶变换程序处理时间抽样信号，并根据实序列离散傅立叶变换的对称性，初步判定程序的正确性；
2.观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT分析信号序列的幅频特性，观察两者的序列形状和频谱曲线有什么异同？
绘出两序列及其幅频特性曲线。
三角波序列
反三角波序列
3.已知余弦信号如下
当信号频率F=50Hz采样间隔T=0.000625s，采样长度N=64时，对该信号进行傅里叶变换。
用FFT程序分析正弦信号，分别在以下情况进行分析。
〔1〕F=50,N=32,T=0.000625；
〔2〕F=50,N=32,T=0.005；
〔3〕F=50,N=32,T=0.0046875；
〔4〕F=50,N=32,T=0.004；
〔5〕F=50,N=64=0.000625
4.选定某一时间信号进行N=64点离散傅立叶变换，并且，对同一信号进行快速傅立叶变换，并比较它们的速度。
三、实验要求
1、编制DFT程序及FFT程序，并比较DFT程序与FFT程序的运行时间。
给编制的程序加注释；
2.完成实验内容2，并对结果进行分析。
在单位圆Z上的变换频谱会相同吗？
如果不同，你能说出那个低频分量更多吗？
为什么？
3.完成实验内容3，并对结果进行分析；
4.利用编制的计算卷积的计算程序，给出一下三组函数的卷积结果。
〔1〕
〔2〕
〔3〕
四、实验程序及其结果分析
1、离散傅里叶变换
function[Xk]=dft（xn,N）%计算N点的DFT
n=0:
1:
N-1;k=0:
1:
N-1;
WN=exp（-j*2*pi/N）;nk=n'*k;
WNnk=WN.^nk;Xk=xn*WNnk;
End
2、〔1〕三角波
clearall;
n=0:
1:
7;
x1=（n+1）.*（n>=0）.*（n<=3）+（8-n）.*（n>=4）;%三角波序列
subplot（3,1,1）;stem（n,x1）;
title（'三角波序列'）;xlabel（'n'）;
N=8;k=0:
N-1;
X1=fft（x1,N）;
magX1=abs（X1）;phaX1=angle（X1）;
subplot（3,1,2）;stem（k,magX1）;
xlabel（'k'）;ylabel（'三角波DFT的幅度'）;
subplot（3,1,3）;stem（k,phaX1）;
xlabel（'k'）;ylabel（'三角波DFT的相位'）;
〔2〕反三角波
clearall;
n=0:
1:
7;
x2=（4-n）.*（n>=0）.*（n<=3）+（n-3）.*（n>=4）;%反三角波序列
subplot（3,1,1）;stem（n,x2）;
title（'反三角波序列'）;xlabel（'n'）;
N=8;k=0:
N-1;
X2=fft（x2,N）;
magX2=abs（X2）;phaX2=angle（X2）;
subplot（3,1,2）;stem（k,magX2）;
xlabel（'k'）;ylabel（'幅度'）;
subplot（3,1,3）;stem（k,phaX2）;
xlabel（'k'）;ylabel（'相位'）;
结果分析：
由图知：
三角波和反三角波序列的波形不同，当N=8时，正反三角波的幅频特性相同，因为两者的时域只差一个相位
3、余弦信号的FFT分析
cleaall;
F=50;N=32;T=0.000625;%
（1）F=50,N=32,T=0.000625
n=1:
N;
x=cos（2*pi*F*n*T）;
figure
（1）;subplot（2,1,1）;plot（n,x）;
ylabel（'x（n）'）;xlabel（'n'）;title（'
（1）F=50,N=32,T=0.000625'）;
X=fft（x）;magX=abs（X）;
subplot（2,1,2）;plot（n,X）;
ylabel（'FFT|X|'）;xlabel（'f（pi）'）;
F=50;N=32;T=0.005;%
（2）F=50,N=32,T=0.005
n=1:
N;
x=cos（2*pi*F*n*T）;
figure
（2）;
subplot（2,1,1）;plot（n,x）;ylabel（'x（n）'）;
xlabel（'n'）;title（'
（2）F=50,N=32,T=0.005'）;
X=fft（x）;magX=abs（X）;
subplot（2,1,2）;plot（n,X）;ylabel（'FFT|X|'）;xlabel（'f（pi）'）;
F=50,N=32,T=0.0046875;%（3）F=50,N=32,T=0.0046875
n=1:
N;
x=cos（2*pi*F*n*T）;
figure（3）;
subplot（2,1,1）;plot（n,x）;ylabel（'x（n）'）;
xlabel（'n'）;title（'（3）F=50,N=32,T=0.0046875'）;
X=fft（x）;
magX=abs（X）;
subplot（2,1,2）;
plot（n,X）;ylabel（'FFT|X|'）;xlabel（'f（pi）'）;
F=50,N=32,T=0.004;%（4）F=50,N=32,T=0.004
n=1:
N;x=cos（2*pi*F*n*T）;
figure（4）;
subplot（2,1,1）;
plot（n,x）;ylabel（'x（n）'）;
xlabel（'n'）;title（'（4）F=50,N=32,T=0.004'）;
X=fft（x）;magX=abs（X）;
subplot（2,1,2）;
plot（n,X）;ylabel（'FFT|X|'）;xlabel（'f（pi）'）;
F=50,N=64,T=0.000625;%（5）F=50,N=64,T=0.000625
n=1:
N;
x=cos（2*pi*F*n*T）;
figure（5）;subplot（2,1,1）;plot（n,x）;ylabel（'x（n）'）;xlabel（'n'）;
title（'（5）F=50,N=64,T=0.000625'）;
X=fft（x）;magX=abs（X）;
subplot（2,1,2）;plot（n,X）;ylabel（'FFT|X|'）;xlabel（'f（pi）'）;
结果分析；不同的采样间隔会产生不同的栅栏效应，相当于透过栅栏欣赏风景,只能看到频谱的一部分,而其它频率点看不见,因此很可能使一部分有用的频率成分被漏掉，从而产生不同的频谱图，减小栅栏效应可用提高采样间隔也就是频率分辨力的方法来解决。
间隔小，频率分辨力高，被“挡住”或丧失的频率成分就会越少。
但会增加采样点数，使计算工作量增加。
4、卷积计算
clearall;
n=0:
14;x1=1.*（n>=0）.*（n<=14）;%第一组函数
h1=（4/5）.^n;
N=27;y=circonvt（x1,h1,N）;
disp（y）;subplot（3,1,1）;
i=0:
N-1;stem（i,y）;ylabel（'式〔1〕卷积'）;
n=0:
9;x2=1.*（n>=0）.*（n<=9）;%第二组函数
n=0:
19;h2=（0.8）.*sin（0.5.*n）;
N=27;y=circonvt（x2,h2,N）;
disp（y）;subplot（3,1,2）;
i=0:
N-1;stem（i,y）;ylabel（'式〔2〕卷积'）;
n=0:
9;x3=（1-0.1.*n）;h3=0.1.*n;%第三组函数
N=17;y=circonvt（x3,h3,N）;
disp（y）;subplot（3,1,3）;
i=0:
N-1;stem（i,y）;ylabel（'式〔3〕卷积'）;
实验三IIR数字滤波器设计
一、试验目的
1、学习模拟－数字变换滤波器的设计方法；
2.掌握双线性变换数字滤波器设计方法；
3.掌握实现数字滤波器的具体方法。
二、实验内容
1、设计一个巴特沃思数字低通滤波器，设计指标如下：
通带内w阻带内w>ws=0.35pi幅度衰减不小于15dB；
2、编制计算设计的数字滤波器幅度特性和相位特性的程序，并进行实验验证。
3、编制实现该数字滤波器程序并且实现数字滤波
〔1〕分别让满足所设计的滤波器的通带、过渡带、阻带频率特性的正弦波通过滤波器，验证滤波器性能；
〔2〕改变正弦抽样时间，验证数字低通滤波器的模拟截止频率实抽样时间的函数。
三、实验要求
1、编制实验内容要求的程序，并给程序加注释；
2、根据实验结果，给出自己设计的数字滤波器的幅度特性和相位特性；
3.用所设计的滤波器对不同频率的正弦波信号进
行滤波，以说明其特性；
4、fp=0.2KHz,Rp=1dB,fs=0.3KHz,As=25dB,T=1ms；分别用脉冲响应不变法及双线性变换法设计一Butterworth数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。
比较这两种方法的优缺点。
四、实验程序及结果分析
1、巴特沃斯低通滤波器
clearall;
wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;
R=10^（-Rp/20）;A=10^（-As/20）;
[b,a]=afd_buttap（wp,ws,Rp,As）;
[C,B,A]=sdir2cas（b,a）;%将直接型装换成级联型
[db,mag,pha,w]=freqs_m（b,a,0.5*pi）;%计算系统频率响应
[ha,x,t]=impulse（b,a）;
subplot（2,1,1）;plot（w/pi,mag）;title（'幅度响应'）;grid;
ylabel（'幅度'）;xlabel（'以\pi为单位的频率'）;
subplot（2,1,2）;plot（w/pi,pha）;title（'相位响应'）;grid;
ylabel（'相位'）;xlabel（'以\pi为单位的频率'）;
2、脉冲响应不变法设计巴特沃斯滤波器
%冲击响应不变法设计数字巴特沃斯低通滤波器
clearall;
T=0.0001;
wp=200*2*pi*T;Rp=1;ws=300*2*pi*T;As=25;%数字滤波器指标
T=1;
omegap=wp/T;omegas=ws/T;%转换为模拟域指标
[cs,ds]=afd_buttap（omegap,omegas,Rp,As）;%模拟巴特沃斯滤波器的计算
%冲击响应不变法
[b,a]=imp_invr（cs,ds,T）;
[C,B,A]=dir2par（b,a）;%将IIR滤波器直接型转换为并联型的函数
%计算数字滤波器频率响应
[db,mag,pha,grd,w]=freqz_m（b,a）;
subplot（3,1,1）;plot（w/pi,mag）;
title（'冲击响应不变法设计巴特沃斯滤波器'）;
axis（[00.501]）;grid;
ylabel（'|幅度|'）;xlabel（'以\pi为单位的频率'）;
subplot（3,1,2）;plot（w/pi,db）;axis（[00.8-4000]）;grid;
ylabel（'对数幅度〔db〕'）;xlabel（'以\pi为单位的频率'）;
3、双线性变换法设计巴特沃斯滤波器
clearall;
T=0.0001;Fs=1/T;
wp=200*2*pi*T;Rp=1;ws=300*2*pi*T;As=25;%数字滤波器指标
%转换成模拟域指标
omegap=（2/T）*tan（wp/2）;
omegas=（2/T）*tan（ws/2）;
%模拟巴特沃斯滤波器的计算
[cs,ds]=afd_buttap（omegap,omegas,Rp,As）;
%双线性变换
[b,a]=bilinear（cs,ds,Fs）;
[C,B,A]=sdir2cas（b,a）;
[db,mag,pha,grd,w]=freqz_m（b,a）;
subplot（3,1,1）;plot（w/pi,mag）;ylabel（'|幅度|'）;title（'双线性变换法设计巴特沃斯滤波器'）;
xlabel（'以\pi为单位的频率'）;axis（[00.501]）;grid;
subplot（3,1,2）;plot（w/pi,db）;axis（[00.8-4000]）;grid;
ylabel（'对数幅度〔db〕'）;xlabel（'以\pi为单位的频率'）;
subplot（3,1,3）;plot（w/pi,pha）;axis（[00.4-44]）;grid;
ylabel（'相位'）;xlabel（'以\pi为单位的频率'）;
4、数字滤波
clearall;
F=0.1;N1=512;T1=0.1;
n=0:
5:
N1;
x=sin（2*pi*F*n*T1）;wp=0.2*pi;ws=0.35*pi;
Rp=1;As=15;T=1;Fs=1/T;
OmegaP=（2/T）*tan（wp/2）;%预畸变
OmegaS=（2/T）*tan（ws/2）;
ep=sqrt（10^（Rp/10）-1）;%求平方根
Ripple=sqrt（1/（1+ep.^2））;Attn=1/10^（As/20）;
N=ceil（（log10（（10^（Rp/10）-1）/（10^（As/10）-1）））/（2*log10（OmegaP/OmegaS）））;
OmegaC=OmegaP/（（10.^（Rp/10）-1）.^（1/（2*N）））;
[cs,ds]=u_buttap（N,OmegaC）;
[b,a]=bilinear（cs,ds,Fs）;
%括号为模拟滤波器系数，中括号中的为数字滤波器系数
[mag,db,pha,w]=freqz_m（b,a）;
y=filter（b,a,x）;
subplot（2,1,1）;stem（x）;
axis（[0100-11]）;title（'输入正弦序列'）;
subplot（2,1,2）;stem（y）;axis（[0100-11]）;
title（'滤波输出序列'）;
结果分析：
滤波输出序列前有一段是由于延时造成的。
实验四FIR数字滤波器
1、FIR滤波器的设计
clearall;
%用哈明窗函数法来设计FIR滤波器选N=51
wp=0.2*pi;ws=0.35*pi;N=51;
n=[0:
1:
N-1];
wc=（ws+wp）/2;%理想低通的截止频率
hd=ideal_lp（wc,N）;%理想低通的冲激响应
w_ham=（hamming（N））';
h=hd.*w_ham%FIR滤波器的冲激响应
[db,mag,pha,grd,w]=freqz_m（h,[1]）;
delta_w=2*pi/1000;
Rp=-（min（db（1:
1:
wp/delta_w+1）））;%实际的通带衰减
As=-round（max（db（ws/delta_w+1:
1:
501）））;%实际的阻带衰减
subplot（2,2,1）;stem（n,hd）;title（'理想冲激响应'）;
axis（[0N-1-0.10.3]）;xlabel（'n'）;ylabel（'hd（n）'）;
subplot（2,2,2）;stem（n,h）;title（'实际冲激响应'）;
axis（[0N-1-0.10.3]）;xlabel（'n'）;ylabel（'h（n）'）;
subplot（2,2,3）;plot（w/pi,pha）;title（'相位响应'）;grid;
xlabel（'以\pi为单位的频率'）;ylabel（'相位'）;
subplot（2,2,4）;plot（w/pi,db）;title（'幅度响应〔db〕'）;grid;
xlabel（'以\pi为单位的频率'）;ylabel（'H（db）'）;
2、分别用矩形窗、汉宁窗、海明窗、三角窗和Blackman窗设计一个线性相位的带通滤波器。
clearall;
w1=0.3*pi;w2=0.5*pi;N=51;
n=[0:
1:
N-1];
hd=ideal_lp（w2,N）-ideal_lp（w1,N）;%理想低通的冲激响应
w_rec=（rectwin（N））';
h1=hd.*w_rec;%矩形窗FIR滤波器的冲激响应
[db,mag,pha,grd,w]=freqz_m（h1,[1]）;
figure;plot（w/pi,db）;title（'矩形窗幅度响应〔db〕'）;
xlabel（'以\pi为单位的频率'）;ylabel（'H1（db）'）;
w_hann=（hann（N））';
h2=hd.*w_hann%汉宁窗FIR滤波器的冲激响应
[db,mag,pha,grd,w]=freqz_m（h2,[1]）;
figure;plot（w/pi,db）;title（'汉宁窗幅度响应〔db〕'）;
xlabel（'以\pi为单位的频率'）;ylabel（'H2（db）'）;
w_ham=（hamming（N））';
h3=hd.*w_ham%哈明窗FIR滤波器的冲激响应
[db,mag,pha,grd,w]=freqz_m（h3,[1]）;
figure;plot（w/pi,db）;title（'海明窗幅度响应〔db〕'）;
xlabel（'以\pi为单位的频率'）;ylabel（'H3（db）'）;
升级会员 