信号处理实验.docx
《信号处理实验.docx》由会员分享，可在线阅读，更多相关《信号处理实验.docx（43页珍藏版）》请在冰豆网上搜索。
信号处理实验
目录
绪论……………………………………………………………………………………………1
1离散时间信号和系统分析
1.1离散时间信号产生与运算…………………………………………………………2
1.2离散时间系统的时域分析…………………………………………………………9
1.3离散时间系统的频域分析………………………………………………………13
1.4离散时间系统频响的零极点确定………………………………………………14
2快速傅立叶变换的应用
2.1FFT的计算………………………………………………………………17
2.2利用FFT进行谱分析…………………………………………………………18
2.3利用FFT实现快速卷积………………………………………………………19
3数字滤波器的设计
3.1数字滤波器的结构…………………………………………………………………23
3.2无限冲激响应（IIR）数字滤波器的设计…………………………………………25
3.3有限冲激响应（FIR）数字滤波器的设计…………………………………………27
4综合应用举例
4.1语音信号处理……………………………………………………………………32
4.2电话拨号音的合成与识别………………………………………………………32
绪论
数字信号处理主要研究如何对信号进行分析、变换、综合、估计与识别等加工处理的基本理论和方法。
随着计算机技术和大规模集成电路技术的发展,数字信号处理以其方便、灵活等特点引起人们越来越多的重视。
在40多年的发展过程中,这门学科基本形成了一套完整的理论体系，其中也包括各种快速、优良的算法,而且数字信号处理的理论和技术也在不断、快速地丰富和完善，新理论和新技术也层出不穷。
学习这门课程的过程中，容易使人感到数字信号处理的概念抽象难懂，其中的分析方法与基本理论不容易很好地理解与掌握。
因此，如何理解与掌握课程中的基本概念、基本原理、基本分析方法以及综合应用所学知识解决实际问题的能力，是本课程学习中所要解决的关键问题。
Matlab是一种面向科学和工程的高级语言，现已成为国际上公认的优秀的科技界应用软件，在世界范围内广为流行和使用。
在欧美高等院校里,Matlab已成为大专院校学生、教师的必要基本技能，广泛应用于科学研究、工程计算、教学等。
上世纪90年代末和本世纪初Matlab在我国也被越来越多地应用于科研和教学工作中。
Matlab是一套功能强大的工程计算及数据处理软件，在工业，电子，医疗和建筑等众多领域均被广泛运用。
它是一种面向对象的，交互式程序设计语言，其结构完整又具有优良的可移植性。
它在矩阵运算，数字信号处理方面有强大的功能。
另外，Matlab提供了方便的绘图功能，便于用户直观地输出处理结果。
本文通过Matlab系列仿真，旨在掌握基本的数字信号处理的理论和方法，提高综合运用所学知识，提高Matlab计算机编程的能力。
进一步加强独立分析问题、解决问题的能力、综合设计及创新能力的培养，同时注意培养实事求是、严肃认真的科学作风和良好的实验习惯。
1.离散时间信号和系统分析
1.1离散时间信号产生与运算
本节的目的是使读者熟悉Matlab中离散时间信号产生和信号运算的基本命令。
几种常用的序列如下：
（1）单位抽样序列
在MATLAB中可以利用zeros（）函数实现：
例如，下列程序
N=input（'Typeinlengthofsequence='）;
n=0:
N-1;
x=zeros（1,N）;
x
（1）=1;
stem（n,x）;
xlabel（'n'）;ylabel（'x（n）'）;
title（'单位抽样序列N取10'）;
输入Typeinlengthofsequence=10，可产生
（2）单位阶越序列
在MATLAB中可以利用ones（）函数实现：
例如，下列程序
N=input（'Typeinlengthofsequence='）;
n=0:
N-1;
x=ones（1,N）;
stem（n,x）;
xlabel（'n'）;ylabel（'x（n）'）;
title（'单位阶越序列N取10'）;
输入Typeinlengthofsequence=10，可产生
（3）正弦序列
在MATLAB中：
例如，下列程序
a=input（'Typeina='）;
b=input（'Typeinb='）;
A=input（'Typeinthegainconstant='）;
N=input（'Typeinlengthofsequence='）;
n=0:
N;x=A*sin（a*pi*n+pi/b）;
stem（n,x）;title（'正弦序列'）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
输入Typeina=0.1，Typeinb=2，Typeinthegainconstant=3，Typeinlengthofsequence=40，可产生
（4）指数序列
在MATLAB中：
例如，下列程序
a=input（'Typeinexponent='）;
K=input（'Typeinthegainconstant='）;
N=input（'Typeinlengthofsequence='）;
n=0:
N;x=K*a.^n;
stem（n,x）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
title（['指数序列alpha=',num2str（a）]）;
输入Typeinexponent=2，Typeinthegainconstant=1，Typeinlengthofsequence=20，可产生如下结果
（5）复指数序列
在MATLAB中：
例如，下列程序
a=input（'Typeinrealexponent='）;
b=input（'Typeinimaginaryexponent='）;
c=a+b*i;
K=input（'Typeinthegainconstant='）;
N=input（'Typeinlengthofsequence='）;
n=1:
N;
x=K*exp（c*n）;
subplot（211）;stem（n,real（x））;
ylabel（'Amplitude'）;
title（'复指数序列Realpart'）;
subplot（212）;stem（n,imag（x））;
xlabel（'Timeindexn'）;
ylabel（'Amplitude'）;
title（'复指数序列Imaginarypart'）;
输入Typeinrealexponent=0.2，Typeinimaginaryexponent=0.2，Typeinthegainconstant=2，Typeinlengthofsequence=40，可产生如下结果
（6）Sinc函数
在MATLAB中：
例如，下列程序
t=-10:
0.01:
10;
x=sinc（t）;
plot（t,x）;
xlabel（'t'）;ylabel（'x（t）'）;
title（'Sinc函数'）;
可产生
（7）随即序列
例如，下列程序
clf;
R=51;
d=0.8*（rand（R,1）-0.5）;
m=0:
R-1;
stem（m,d','b'）;
title（'随机序列'）;
xlabel（'k'）;ylabel（'f（k）'）;
可产生
序列的基本运算有：
（1）序列加法和乘法
在MATLAB中：
x=c+b；
y=c.*b；
例如，下列程序
%取a=[2，1,3,4]，b=[0，1，2,3,1]
m=1:
4;
a=[2134];
c=[21340];
n=1:
5;
b=[01231];
c=[azeros
（1）];
x=c+b;
y=c.*b;
subplot（4,1,1）;stem（m,a）;xlabel（'m'）;ylabel（'a（m）'）;
subplot（4,1,2）;stem（n,b）;xlabel（'n'）;ylabel（'b（n）'）;
subplot（4,1,3）;stem（n,x）;xlabel（'n'）;ylabel（'x（n）'）;
title（'序列的加法'）;
subplot（4,1,4）;stem（n,y）;xlabel（'n'）;ylabel（'y（n）'）;
title（'序列的乘法'）;
可产生
（2）序列的卷积
在MATLAB中：
c=conv（a,b）；
例如，下列程序
a=input（'Typeinthefirstsequence='）;
b=input（'Typeinthesecondsequence='）;
c=conv（a,b）;
M=length（c）-1;n=0:
1:
M;
disp（'outputsequence='）;disp（c）；stem（n,c）；
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;title（'序列的卷积'）;
输入Typeinthefirstsequence=[123]，Typeinthesecondsequence=[456]，可产生：
outputsequence=
413282718
1.2离散时间系统的时域分析
对线性离散时间系统，若y1[n]和y2[n]分别是输入序列x1[n]和x2[n]的响应，则输入
x[n]=ax1[n]+bx2[n]
的输出响应为
y[n]=ay1[n]+by2[n]
式中叠加性质对任意常数a和b以及任意输入x1[n]和x2[n]都成立。
反之，则系统称之为非线性。
例如，下列程序
%y[n]-0.4y[n-1]+0.75y[n-2]=2.2403x[n]+2.4908x[n-1]+2.2403x[n-2]
n=0:
40;
a=2;b=-3;
x1=cos（2*pi*0.1*n）;
x2=sin（2*pi*0.1*n）;
x=a*x1+b*x2;
num=[2.24032.49082.2403];den=[1-0.40.75];
ic=[00];%设置零初始条件
y1=filter（num,den,x1,ic）;%计算输出y1[n]
y2=filter（num,den,x2,ic）;%计算输出y2[n]
y=filter（num,den,x,ic）;%计算输出y[n]
yt=a*y1+b*y2;
d=y-yt;%计算差值输出d[n]
%画出输出和差信号
subplot（3,1,1）;stem（n,y）;ylabel（'振幅'）;
title（'加权输入:
a\cdotx_{1}[n]+b\cdotx_{2}[n]的输出'）;
subplot（3,1,2）;stem（n,yt）;ylabel（'振幅'）;
title（'加权输出:
a\cdoty_{1}[n]+b\cdoty_{2}[n]'）;
subplot（3,1,3）;stem（n,d）;
xlabel（'时间序号n'）;ylabel（'振幅'）;title（'差信号'）;
可产生
对于离散时不变系统，若y1[n]是x1[n]的响应，则输入
x[n]=x1[n-n0]
的输出响应为
y[n]=y1[n-n0]
式中n0时任意整数。
上面的输入输出关系，对任意输入序列及其相应的输出成立。
反之，则系统称之为时变的。
例如，下列程序
%y[n]-0.4y[n-1]+0.75y[n-2]=2.2403x[n]+2.4908x[n-1]+2.2403x[n-2]
clf;
n=0:
40;D=10;a=3.0;b=-2;
x=a*cos（2*pi*0.1*n）+b*sin（2*pi*0.1*n）;
xd=[zeros（1,D）x];
num=[2.24032.49082.2403];
den=[1-0.40.75];
ic=[00];%设置零初始条件
y=filter（num,den,x,ic）;%计算输出y[n]
yd=filter（num,den,xd,ic）;%计算输出yd[n]
d=y-yd（1+D:
41+D）;%计算差值输出d[n]
%画出输出
subplot（3,1,1）;
stem（n,y）;ylabel（'振幅'）;
title（'输出y[n]'）;grid;
subplot（3,1,2）;
stem（n,yt（1:
41））;ylabel（'振幅'）;
title（['由于延时输入x[n',num2str（D）,']的输出']）;grid;
subplot（3,1,3）;
stem（n,d）;
xlabel（'时间序号n'）;ylabel（'振幅'）;
title（'差值信号'）;grid;
可产生结果
离散时间系统的仿真：
线性和非线性系统、时变和非时变系统的仿真离散系统
其输入、输出关系可用以下差分方程描述：
输入信号分解为冲激信号
记系统单位冲激响应
，则系统响应为如下的卷积计算式：
当
时，h[n]是有限长度的（n：
[0，M]），称系统为FIR系统；反之，称系统为IIR系统。
1.3离散时间系统的频域分析
序列x[n]的DTFT定义：
在MATLAB中，可用freqz计算出离散时间系统的频率响应。
可用下列程序计算差分方程
y（n）+0.7y（n-1）-0.45y（n-2）-0.6y（n-3）=0.8x（n）-0.44x（n-1）+0.36x（n-2）+0.02x（n-3）
的单位脉冲响应：
%x（n）=zeros（1,N-1）,0<=n<=40
N=41;
a=[0.8-0.440.360.22];b=[10.7-0.45-0.6];
x=[1zeros（1,N-1）];
k=0:
1:
N-1;y=filter（a,b,x）;
stem（k,y）;
xlabel（'n'）;ylabel（'幅度'）;
可产生
可用下列程序计算差分方程
y（n）+0.7y（n-1）-0.45y（n-2）-0.6y（n-3）=0.8x（n）-0.44x（n-1）+0.36x（n-2）+0.02x（n-3）
对应系统函数的DTFT：
k=256;
num=[0.8-0.440.360.02];
den=[10.7-0.45-0.6];
w=0:
pi/k:
pi;
h=freqz（num,den,w）;
subplot（2,2,1）;plot（w/pi,real（h））;grid;
title（'实部'）;xlabel（'\omega/\pi'）;ylabel（'幅度'）;
subplot（2,2,2）;plot（w/pi,imag（h））;grid;
title（'虚部'）;xlabel（'\omega/\pi'）;ylabel（'幅值'）;
subplot（2,2,3）;plot（w/pi,abs（h））;grid;
title（'幅度谱'）;xlabel（'\omega/\pi'）;ylabel（'幅值'）;
subplot（2,2,4）;plot（w/pi,angle（h））;grid;
title（'相位谱'）;xlabel（'\omega/\pi'）;ylabel（'弧度'）;
可产生
1.4离散时间系统频响的零极点确定
离散系统的时域方程为
其变换域分析方法如下：
（1）频域
系统的频率响应为
（2）Z域
系统的转移函数为
分解因式
其中
和
称为零、极点。
在MATLAB中，可以用函数[z，p，K]=tf2zp（num，den）求得有理分式形式的系统转移函数的零、极点，用函数zplane（z，p）绘出零、极点分布图；也可以用函数zplane（num，den）直接绘出有理分式形式的系统转移函数的零、极点分布图。
可用下列程序，求解已知离散系统H（z）的零极点图，并求解h（k）和H（e^jw）:
b=[121];
a=[1-0.5-0.0050.3];
subplot（311）;
zplane（b,a）;
axis（[-3,3,-1,1]）
num=[0121];
den=[1-0.5-0.0050.3];
h=impz（num,den）;
subplot（312）;
stem（h）;
xlabel（'k'）;
ylabel（'h（k）'）;
[H,w]=freqz（num,den）;
subplot（313）;
plot（w/pi,abs（H））;
xlabel（'/omege'）;
ylabel（'abs（H）'）;
可产生系统H（z）的零极点图，以及h（k）和H（e^jw）：
2FFT的应用
2.1快速傅立叶变换的计算
N点序列的DFT和IDFT变换定义式如下：
利用旋转因子
具有周期性，可以得到快速算法（FFT）。
在MATLAB中，可以用函数U=fft（u，N）和u=ifft（U，N）计算N点序列的DFT正、反变换。
例如，可用下列程序求x=cos（5*pi*n/16）的16点序列的16和32点的DFT
N1=16;N2=32;
n1=0:
N1-1;n2=0:
N2-1;
a=cos（5*pi*n1/16）;b=cos（5*pi*n2/16）;
x1=fft（a,N1）;x2=fft（b,N2）;
subplot（2,1,1）;stem（n1,abs（x1）,'.'）;axis（[0,20,0,20]）;
xlabel（'n1'）;ylabel（'|X1（n1）|'）;title（'16点DFT'）;
subplot（2,1,2）;stem（n2,abs（x2）,'.'）;axis（[0,20,0,20]）;
xlabel（'n2'）;ylabel（'|X2（n2）|'）;title（'32点DFT'）;
可得到
2.2利用FFT进行谱分析
对信号进行谱分析，就是计算信号的傅里叶变换。
用Matlab语言编制信号产生子程序，产生典型信号供谱分析用，其中x（n）是由两个正弦信号及白噪声的叠加，产生两个正弦加白噪声，对产生的信号进行谱分析，绘出序列和幅频特性曲线。
对连续信号进行谱分析也是连续的，应先对连续信号进行时域采样，再利用DFT对采样序列进行谱分析。
对信号进行谱分析，就是计算信号的傅里叶变换，绘出序列和幅频特性曲线。
然而，因为对连续信号进行谱分析也是连续的，应先对连续信号进行时域采样，故而会使谱分析引入误差，所以用DFT对连续信号进行谱分析的结果都是近似的。
例如，下列程序是利用FFT对
x（n）=a1*sin（w*f1（0:
N-1））+sin（w*f2*（0:
N-1））+randn（1,N）w=2*pi/fs
进行谱分析,x（n）是由两个正弦信号及白噪声的叠加
N=256;a1=5;a2=3;
f1=.1;f2=.2;fs=1;
w=2*pi/fs;
x=a1*sin（w*f1*（0:
N-1））+sin（w*f2*（0:
N-1））+randn（1,N）;
%应用FFT求频谱
subplot（2,1,1）;plot（x（1:
N/4））;title（'原始信号'）;
y=fft（x）;
subplot（2,1,2）;plot（f,（abs（y）））;title（'频域信号'）;
可得到
2.3利用FFT实现快速卷积
运用DFT的快速算法FFT，对序列进行卷积，当N很大时，计算速度会快很多，用FFT计算线形卷积步骤为：
（1）求H（k）=FFT[h（n）],N点；
（2）求X（k）=DFT[x（n）],N点；
（3）计算Y（k）=X（k）H（k）；
（4）求y（n）=IDFT[Y（k）],N点。
下列程序可实现序列xn=0.8.^n和hn=ones（1,N2）的快速卷积，
n=[0:
1:
50];
m=[0:
1:
20];
N1=length（n）;
N2=length（m）;
xn=0.8.^n;%生成x（n）
hn=ones（1,N2）;%生成h（n）
N=N1+N2-1;
XK=fft（xn,N）;
HK=fft（hn,N）;
YK=XK.*HK;
yn=ifft（YK,N）;
ifall（imag（xn）==0）&（all（imag（hn）==0））%实序列的循环卷积仍然为实序列
yn=real（yn）;
end
x=0:
N-1;
stem（x,yn,'.'）;
得到卷积序列：
3数字滤波器的设计
3.1数字滤波器的结构
无限冲击响应传输函数为
有限冲激响应传输函数为
在MATLAB中无限冲激响应传输函数可用[z,p,k]=tf2zp（num,den）实现，有限冲激响应传输函数可用[r1,p1,k1]=residuez（num,den）和[r2,p2,k2]=residue（num,den）实现。
下列程序，可实现一个无限冲激响应传输函数的并联形式
num=input（'分子系数向量='）;
den=input（'分母系数向量='）;
[r1,p1,k1]=residuez（num,den）;
[r2,p2,k2]=residue（num,den）;
disp（'并联I型'）;
disp（'零点是'）;disp（r1）;
disp（'极点是'）;disp（p1）;
disp（'常数'）;disp（k1）;
disp（'并联II型'）;
disp（'零点是'）;disp（r2）;
disp（'极点是'）;disp（p2）;
disp（'常数'）;disp（k2）;
运行后输入：
分子系数向量=[123]，分母系数向量=[456]，可得
并联I型
零点是
-0.1250-0.0148i
-0.1250+0.0148i
极点是
-0.6250+1.0533i
-0.6250-1.0533i
常数
0.5000
并联II型
零点是
0.0938-0.1224i
0.0938+0.1224i
极点是
-0.6250+1.0533i
-0.6250-1.0533i
常
升级会员 