数字信号处理毕业课程设计实验报告.docx
《数字信号处理毕业课程设计实验报告.docx》由会员分享，可在线阅读，更多相关《数字信号处理毕业课程设计实验报告.docx（35页珍藏版）》请在冰点文库上搜索。
数字信号处理毕业课程设计实验报告
（此文档为word格式，下载后您可任意编辑修改！
）
基础实验
实验一离散时间系统及离散卷积
一、实验原理
利用Matlab软件计算出系统函数的零极点分布、单位脉冲响应和系统频率响应等的图像并于笔算结果进行比较，找出异同。
编译合适程序能计算取值范围不同的离散卷积。
二、实验目的
（1）熟悉MATLAB软件的使用方法。
（2）熟悉系统函数的零极点分布、单位脉冲响应和系统频率响应等概念。
（3）利用MATLAB绘制系统函数的零极点分布图、系统频率响应和单位脉冲响应。
（4）熟悉离散卷积的概念，并利用MATLAB计算离散卷积。
三、实验源程序及实验结果
1.离散时间系统的单位脉冲响应
functionpr1（）%定义函数pr1
a=[1,-1,0.9];%定义差分方程y（n）-y（n-1）+0.9y（n-2）=x（n）
b=1;
x=impseq（）;%调用impseq函数
n=[-20:
120];%定义n从-20到120
ylabel（'h（n）'）;%绘图纵座标为h（n）
figure
（2）%绘图figure2
[z,p,g]=tf2zp（b,a）;%绘出零极点图
zplane（z,p）
function[x,n]=impseq（n0,n1,n2）%声明impseq函数
n=[n1:
n2];
x=[（n-n0）==0];
2.离散系统的幅频、相频的分析方法
functionpr2（）
b=[0.0181,0.0543,0.0543,0.0181];
a=[1.000,-1.76,1..2781];
m=0:
length（b）-1;%m从0到3
l=0:
length（a）-1;%l从0到3
K=500;k=1:
K;w=pi*kK;%角频率w
H=（b*exp（-j*m'*w））.（a*exp（-j*l'*w））;%对系统函数的定义
magH=abs（H）;%magH为幅度
angH=angle（H）;%angH为相位
figure
（1）
subplot（2,1,1）;%在同一窗口的上半部分绘图
plot（wpi,magH）;%绘制w（pi）-magH的图形
grid;
axis（[]）;%限制横纵座标从0到1
xlabel（'w（pi）'）;%x座标为w（pi）
ylabel（'|H|'）;%y座标为angle（H）
title（'幅度，相位响应'）;%图的标题为:
'幅度，相位响应'
subplot（2,1,2）;%在同一窗口的下半部分绘图
plot（wpi,angH）;%绘制w（pi）-angH的图形
grid;%为座标添加名称
xlabel（'w（pi）'）;%x座标为w（pi）
ylabel（'angle（H）'）;%y座标为angle（H）
3.离散卷积的计算
functionpr3（）
n=-5:
50;%声明n从-5到50
u1=stepseq（）;%调用stepseq函数声用明u1=u（n）
u2=stepseq（）;%调用stepseq函数声用明u2=u（n-10）
x=u1-u2;%x（n）=u（n）-u（n-10）
）.*u1;%）=0.9^n*u（n）
L=length（x）+length（,x）;axis（[]）;title（'输入序列'）;xlabel（'n'）;ylabel（'x（n）'）;%输入序列
subplot（3,1,2）;stem（n,）'）;%冲激相应序列
n=0:
L-1;
subplot（3,1,3）;stem（n,real（y））;title（'输出响应'）;xlabel（'n'）;ylabel（'y（n）'）;%卷积结果
function[x,n]=stepseq（n0,n1,n2）
n=n1:
n2;
x=[（n-n0）>=0];
实验二离散傅立叶变换与快速傅立叶变换
一、实验原理
对有限长序列使用离散Fouier变换（DFT）可以很好的反映序列的频谱特性，而且易于用快速算法在计算机上实现，当序列x（n）的长度为N时，它的DFT定义为
反变换为
有限长序列的DFT是其Z变换在单位圆上的等距采样，或者说是序列Fourier变换的等距采样，因此可以用于序列的谱分析。
FFT是为了减少DFT运算次数的一种快速算法。
它是对变换式进行一次次分解，使其成为若干较短序列的组合，从而减少运算量。
常用的FFT是以2为基数的，其长度。
它的效率高，程序简单，使用非常方便，当要变换的序列长度不等于2的整数次方时，为了使用以2为基数的FFT，可以用末位补零的方法，使其长度延长至2的整数次方。
用FFT可以实现两个序列的圆周卷积。
在一定的条件下，可以使圆周卷积等于线性卷积。
一般情况，设两个序列的长度分别为N1和N2，要使圆周卷积等于线性卷积的充要条件是FFT的长度
N≥N1＋N2
对于长度不足N的两个序列，分别将他们补零延长到N。
二、实验目的
（1）加深理解离散傅立叶变换及快速傅立叶变换概念；
（2）学会应用FFT对典型信号进行频谱分析的方法；
（3）研究如何利用FFT程序分析确定性时间连续信号；
（4）熟悉应用FFT实现两个序列的线性卷积的方法。
3、实验源程序及结果
1.DFT程序
functionpr4（）
F=50;N=64;T=0.000625;n=1:
N;
x=cos（2*pi*F*n*T）;%x（n）=cos（pi*n16）
subplot（2,1,1）;stem（n,x）;title（'x（n）'）;xlabel（'n'）;%x（n）
X=dft（x,N）;
subplot（2,1,2）;stem（n,X）;title（'DFT|X|'）;xlabel（'f（pi）'）;%DFT|X|
%dft的子程序
%实现离散傅里叶变换
function[Xk]=dft（xn,N）
n=0:
N-1;
k=0:
N-1;
WN=exp（-j*2*piN）;
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk;
2.三角序列与反三角序列DFT
N=8;
fs=100;T=1fs;
n=-6:
9;
xc=（n+1）.*（n>=0&n<=3）+（8-n）.*（n>=4）;
xd=（4-n）.*（n>=0&n<=3）+（n-3）.*（n>=4）;
X1=fft（xc,N）;
X2=fft（xd,N）;
magX1=abs（X1）,phaX1=angle（X1）
magX2=abs（X2）,phaX2=angle（X2）
subplot（3,2,1）;stem（n,xc）;
xlabel（'n'）;title（'xc（n）'）;
subplot（3,2,2）;stem（n,xd）;
xlabel（'n'）;title（'xd（n）'）;
k=0:
N-1;
subplot（3,2,3）;stem（k,magX1）;
subplot（3,2,4）;stem（k,magX2）;
xlabel（'k'）;ylabel（'dftdefudu'）;
subplot（3,2,5）;stem（k,phaX1）;
subplot（3,2,6）;stem（k,phaX2）;
xlabel（'k'）;ylabel（'dftdexiangwei'）;
3.余弦信号傅里叶变换
F=50;N=64;T=0.000625;n=1:
N;
x=cos（2*pi*F*n*T）;%x（n）=cos（pi*n16）
subplot（2,1,1）;plot（n,x）;title（'x（n）'）;xlabel（'n'）;%在第一个子窗中绘图x（n）
X=fft（x）;
subplot（2,1,2）;plot（n,X）;title（'DTFT|X|'）;xlabel（'f（pi）'）;%在第二个子图中绘图x（n）的快速傅里叶变换
4.计算卷积
functionpr6（）
n=0:
15;x=1.^n;;
x（16:
32）=0;%x（n）=1,n=0~15;x（n）=0,n=16~32
）=（0.8）^n,n=0~15;）=0,n=16~32
subplot（3,1,1）;stem（x）;title（'x（n）'）;axis（[.5]）;%x（n）
subplot（3,1,2）;stem（）'）;axis（[.5]）;%）
X=fft（x）;%X（n）为x（n）的快速傅里叶变换
H=fft（）为h（n）的快速傅里叶变换
Y=X.*H;%Y（n）=X（n）*H（n）
y=ifft（Y）;%y（n）为Y（n）的傅里叶反变换
subplot（3,1,3）;stem（abs（y））;title（'y（n=x（n）*））'）;axis（[]）;%y（n）
四、实验总结与思考
1、在较短的傅里叶变换中，FFT的计算速度与DFT相比不是很明显，序列计算长度越长，计算时间差距越大，FFT较快；
2、对于不同序列的较小长度的频谱分析可能会得到相同的频谱，适当加倍长度会避免这种情况的发生；
3、对同一序列的不同间隔的FFT变换，在满足奈奎斯特定律的情况下也会产生栅栏效应、频谱泄露、旁瓣效应等，采取适当的方法可以减弱这些不利效应；
4、在计算两个序列的离散卷积的时候要注意序列的长度L>=M+N-1。
实验三IIR数字滤波器设计
一、实验原理
（1）脉冲响应不变法
用数字滤波器的单位脉冲响应序列模仿模拟滤波器的冲激响应,让正好等于的采样值，即其中T为采样间隔，如果以及分别表示的拉氏变换及的Z变换，则
（2）双线性变换法
s平面与z平面之间满足以下映射关系：
以低通数字滤波器为例，将设计步骤归纳如下：
1.确定数字滤波器的性能指标：
通带临界频率fp、阻带临界频率fs；通带内的最大衰减Rp；阻带内的最小衰减As；采样周期T；
2.确定相应的数字角频率，ωp=2πfpT；ωs=2πfrT；
3.计算经过预畸的相应模拟低通原型的频率，
根据Ωp和Ωs计算模拟低通原型滤波器的阶数N，并求得低通原型的传递函数Ha（s）；
4.用上面的双线性变换公式代入Ha（s），求出所设计的传递函数H（z）；分析滤波器特性，检查其指标是否满足要求。
二、实验目的
（1）学习模拟－数字变换滤波器的设计方法；
（2）掌握双线性变换数字滤波器设计方法；
（3）掌握实现数字滤波器的具体方法。
三、实验源程序及结果
functionpr7（）
wp=0.2*pi;ws=0.3*pi;
Rp=1;As=15;T=1;Fs=1T;
OmegaP=（2T）*tan（wp2）;%OmegaP（w）=2*tan（0.1*pi）
OmegaS=（2T）*tan（ws2）;%OmegaS（w）=2*tan（0.15*pi）
ep=sqrt（10^（Rp10）-1）;
Ripple=sqrt（1（1+ep.^2））;
Attn=110^（As20）;
N=ceil（（log10（（10^（Rp10）-1）（10^（As10）-1）））（2*log10（OmegaPOmegaS）））;
OmegaC=OmegaP（（10.^（Rp10）-1）.^（1（2*N）））;
[cs,ds]=u_buttap（N,OmegaC）;
[b,a]=bilinear（cs,ds,Fs）;
[mag,db,pha,w]=freqz_m（b,a）;
%在第一个子窗绘制幅度响应的图形
subplot（3,1,1）;plot（wpi,mag）;title（'幅度响应'）;xlabel（'w（pi）'）;ylabel（'H'）;
axis（[.1]）;
set（gca,'XTickmode','manual','XTick',[0,0.2,0.3,1.1]）;
set（gca,'YTickmode','manual','YTick',[0,Attn,Ripple,1]）;
grid;
%在第二个子窗以分贝为单位绘制幅度响应的图形
subplot（3,1,2）;plot（wpi,db）;title（'幅度响应（dB）'）;xlabel（'w（pi）'）;ylabel（'H'）;
axis（[]）;
set（gca,'XTickmode','manual','XTick',[0,0.2,0.3,1.1]）;
set（gca,'YTickmode','manual','YTick',[]）;
grid;
%在第三个子窗绘制相位响应的图形
subplot（3,1,3）;plot（wpi,pha）;title（'相位响应'）;xlabel（'w（pi）'）;ylabel（'piunit'）;
set（gca,'XTickmode','manual','XTick',[0,0.2,0.3,1.1]）;
set（gca,'YTickmode','manual','YTick',[-1,0,1]）;
grid;
%非归一化巴特沃斯模拟低通滤波器原型函数
function[b,a]=u_buttap（N,OmegaC）
[z,p,k]=buttap（N）;
p=p*OmegaC;k=k*OmegaC.^N;B=real（poly（z））;
b0=k;b=k*B;a=real（poly（p））;
%计算模拟滤波器频率响应的函数
function[mag,db,pha,w]=freqz_m（b,a）
[H,w]=freqz（b,a,1000,'whole'）;
H=（H（1:
501））';w=（w（1:
501））';
mag=abs（H）;pha=angle（H）;
db=20*log10（（mag+eps）max（mag））;
四、实验总结与思考
1、数字滤波器可以用模拟滤波器通过模拟到数字的变换来实现，具体方法有三种：
冲激响应不变法、阶跃响应不变法和双线性变换法等；
2、冲激响应不变法仅适合于基本上是限带的低通滤波器，对于高通或带阻滤波器应该附加限带要求，以避免严重的混叠失真；
3、双线性变换法可以用频率预畸变的方法来补偿频率标度的非线性失真；
4、数字低通滤波器的模拟截止频率是抽样时间的函数；
5、低通滤波器对相应频率的输入信号有相应限制作用。
实验四FIR数字滤波器设计
一、实验原理
窗函数法设计线性相位FIR滤波器步骤：
∙确定数字滤波器的性能要求：
临界频率{ωk}，滤波器单位脉冲响应长度N；
∙根据性能要求，合理选择单位脉冲响应的奇偶对称性，从而确定理想频率响应的幅频特性和相频特性；
∙求理想单位脉冲响应，在实际计算中，可对按M（M远大于N）点等距离采样，并对其求IDFT得，用代替；
∙选择适当的窗函数，根据求所需设计的FIR滤波器单位脉冲响应；
∙求，分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。
窗函数的傅式变换的主瓣决定了过渡带宽。
的旁瓣大小和多少决定了在通带和阻带范围内波动幅度，常用的几种窗函数有：
∙矩形窗w（n）=RN（n）
∙Hanning窗
∙Hamming窗
∙Blackmen窗
二、实验目的
（1）学习FIR数字滤波器窗口函数设计法；
（2）熟悉线性相位FIR滤波器的幅频特性和相频特性；
（3）了解各种不同窗函数对滤波器性能的影响；
（4）进行FIR、IIR滤波器的性能比较。
3、实验源程序及运行结果
1.低通滤波器相位响应‘幅度响应
functionpr8（）
wp=0.2*pi;ws=0.3*pi;
tr_width=ws-wp;
M=51;
n=0:
M-1;%窗口长度n
wc=（ws+wp）2;%理想低通的截止频率
alpha=（M-1）2;
m=n-alpha+eps;%加一个很小的数以避免用零来除
（wc*m）.（pi*m）;%理想低通的冲激响应
w_（db（1:
wpdelta_w+1）））;%实际的通带衰减
As=-round（max（db（wsdelta_w+1:
501）））;%实际的最小阻带衰减
subplot（2,2,1）;stem（n,）'）;
subplot（2,2,2）;stem（n,）'）;
subplot（2,2,3）;plot（wpi,pha）;title（'滤波器相位响应'）;
axis（[0,1,-pi,pi]）;ylabel（'pha'）;
set（gca,'XTickmode','manual','XTick',[0,0.2,0.3,1.1]）;
set（gca,'YTickmode','manual','YTick',[-pi,0,pi]）;
grid;
subplot（2,2,4）;plot（wpi,db）;title（'滤波器幅度响应'）;
axis（[]）;ylabel（'H（db）'）;
set（gca,'XTickmode','manual','XTick',[0,0.2,0.3,1.1]）;
set（gca,'YTickmode','manual','YTick',[]）;
function[mag,db,pha,w]=freqz_m（b,a）
[H,w]=freqz（b,a,1000,'whole'）;
H=（H（1:
501））';w=（w（1:
501））';
mag=abs（H）;pha=angle（H）;
db=20*log10（（mag+eps）max（mag））;
2.分别用矩形窗、汉宁窗、海明窗、三角窗和Blackman窗设计线性相位带通滤波器
functionpr8（）
wp=0.3*pi;
ws=0.5*pi;
tr_width=ws-wp;
M=ceil（6.6*pitr_width）+1;
n=0:
M-1;
wc=（ws+wp）2;
alpha=（M-1）2;
m=n-alpha+eps;
（ws*m）.（pi*m）-sin（wp*m）.（pi*m）;
w_（M））';
（db（1:
wpdelta_w+1）））;
As=-round（max（db（wsdelta_w+1:
501）））;
subplot（2,2,1）;
stem（n,）'）;
subplot（2,2,2）;
stem（n,）'）;
subplot（2,2,3）;
plot（wpi,pha）;
title（'滤波器相位响应'）;
axis（[0,1,-pi,pi]）;
ylabel（'pha'）;
set（gca,'XTickmode','manual','XTick',[0,0.2,0.3,1.1]）;
set（gca,'YTickmode','manual','YTick',[-pi,0,pi]）;
grid;
subplot（2,2,4）;
plot（wpi,db）;
title（'滤波器幅度响应'）;
axis（[]）;
ylabel（'H（db）'）;
set（gca,'XTickmode','manual','XTick',[0,0.2,0.3,1.1]）;
set（gca,'YTickmode','manual','YTick',[]）;
function[mag,db,pha,w]=freqz_m（b,a）
[H,w]=freqz（b,a,1000,'whole'）;
H=（H（1:
501））';
w=（w（1:
501））';
mag=abs（H）;
db=20*log10（（mag+eps）max（mag））;
pha=angle（H）;
矩形窗
汉宁窗
哈明窗
三角窗
Blackman窗
四、实验总结与思考
1、IIR滤波器的系统函数是有理分式，其分母多项式对应于反馈支路，因而这种滤波器是递归结构的系统，只有当所有极点都在单位圆内时滤波器才是稳定的；
2、IIR滤波器可以用比FIR滤波器少的阶数来满足相同的技术指标，这样，IIR滤波器所用的存储单元和所需的运算次数都比FIR滤波器少；
3、IIR滤波器可利用模拟滤波器现成的设计公式、数据和表格，计算工作量较小，对计算工具要求不高；FIR滤波器没有现成的设计公式，窗函数法只给出窗函数的计算公式，但计算通带和阻带衰减仍无显式表达式，一般，FIR滤波器的设计只有计算机程序可以利用；
4、FIR滤波器的相位严格线性；
5、相同指标情况下，FIR滤波器的系统特性更好些；
6、改变FIR滤波器的阶数，对相位特性的影响比幅度特性更显著；
7、如果需要严格线性相位的滤波器，则需要用FIR滤波器。
综合实验
一、实验目的
借助基础实验篇编制的程序，对语音信号进行处理。
通过综合实验篇，充分了解信号处理及数字信号处理的过程；使得学生能够对自己设计的频谱分析程序与滤波器程序得到灵活的应用；使得学生能够综合分析信号处理过程各个阶段信号时域与频域特性。
二、实验内容
1、录制一段各人自己的语音信号。
2、对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；
3、给定滤波器的性能指标，采用窗函数法或双线性变换设计滤波器，并画出滤波器的频率响应。
三、实验源程序及运行结果
fs=20000;
[x]=wavread（'E:
\张利飞'）;
n=length（x）;
Noise=0.2*randn（n,2）;
x2=x+Noise;
%sound（x）;
Ts=1fs;R1=10;
wp=2*pi*1000fs;
ws=2*pi*1200fs;
Rp=1;
Rl=100;
wp1=2Ts*tan（wp2）;%将模拟指标转换成数字指标
ws1=2Ts*tan（ws2）;
[N,Wn]=buttord（wp1,ws1,Rp,R1,'s'）;%选择滤波器的最小阶数
[Z,P,K]=buttap（N）;%创建butterworth模拟滤波器
[Bap,Aap]=zp2tf（Z,P,K）;
[b,a]=lp2lp（Bap,Aap,Wn）;
[bz,az]=bilinear（b,a,fs）;%用双线性变换法实现模拟滤波器到数字滤波器的转换
[H,W]=freqz（bz,az）;%绘制频率响应曲线
figure
（1）
plot（W*fs（2*pi）,abs（H））;
title（'IIR低通滤波器'）;
grid;
axis（[0300002]）;
xlabel（'频率Hz'）
ylabel（'频率响应幅度'）
f1=filter（bz,az,x2）;
figure
（2）
subplot（2,1,1）
plot（x2）%画出滤波前的时域图
title（'IIR低通滤波器滤波前的时域波形'）;
subplot（2,1,2）
plot（f1）;%画出滤波后的时域图
title（'IIR低通滤波器滤波后的时域波形'）;
sound（f）;%播放滤波后的信号
F0=fft（f1,1024）;
f=fs*（0:
511）1024;
figure（3）
y2=fft（x2,1024）;
subplot（2,1,1）;
升级会员 