离散信号和系统分析的MATLAB实现.docx
《离散信号和系统分析的MATLAB实现.docx》由会员分享，可在线阅读，更多相关《离散信号和系统分析的MATLAB实现.docx（50页珍藏版）》请在冰豆网上搜索。
离散信号和系统分析的MATLAB实现
1.9
离散信号和系统分析的MATLAB实现
1.9.1利用MATLAB产生离散信号
用MATLAB表示一离散序列x[k]时，可用两个向量来表示。
其中一个向量表示自变量k的取值范围，另一个向量表示序列x[k]的值。
例如序列x[k]={2,1,1,-1,3,0,2}可用MATLAB表示为
K=-2:
4;x=[2,1,1,-1,3,0,2]
可用stem（k,f）画出序列波形。
当序列是从k=0开始时，可以只用一个向量x来表示序列。
由于计算机内寸的限制，MATLAB无法表示一个无穷长的序列。
例1-38利用MATLAB计算单位脉冲序
范围内各点的取值。
解：
%progran1_1产生单位脉冲序列
Ks=-4;ke=4;n=2;
K=[ks:
ke];
X=[（k-n）==0];
Stem（k,x）:
xlabel（‘k’）;
程序产生的序列波形如图1-49所示。
例1-39利用MATLAB画出信号
X
）+n[k],
的波形。
其中n[k]表示为均值为0方差为1的Gauss分布随机信号。
解：
MALAB提供了两个产生（伪）随机序列的函数。
Rand（1,N）产生1行N列的[0，1]均匀分布随机数。
Randn（1,N）产生1行N列均值为0方差为1的Gauss分布随机数。
%program1_2产生受噪声干扰的正弦信号
N=100;k=0:
N;
X=10*sin（0.02*pi*k）+randn（1,N+1）;
Plot（k,x）;
Xlabel（‘k’）;
Ylabel（‘x[k]’）;
程序产生序列如图1-50所示。
1.9.2离散卷积的计算
离散卷积是数字信号处理中的一个基本运算，MTLAB提供的计算两个离散序列卷积的函数是conv，其调用方式为
y=conv（x,h）
其中调用参数x,h为卷积运算所需的两个序列，返回值y是卷积结果。
MATLAB函数conv的返回值y中只有卷积的结果，没有y的取值范围。
由离散序列卷积的性质可知，当序列x和h的起始点都为k=0时，y的取值范围为k=0至length（x）+length（h）-2。
当序列x或（和）h的起始点不是k=0时，由例1-3知，y的非零值范围可根据例1-4的结论进行计算。
例1-40试用MATLAB函数conv计算例1-2中序列的卷积。
解：
program1_3计算离散卷积
x=[-0.5,0,0.5,1];%序列x的值
kx=-1:
2;%序列x的取值范围
h=[1,1,1];
kh=-2:
0;
y=conv（x,h）;%计算卷积
k=kx
（1）+kh
（1）:
kx（end）+kh（end）;%计算y的取值范围
stem（k,y）;
xlabel（‘k’）;ylabel（‘y’）;
程序的运行结果如图1-51所示。
离散LTI系统响应MATLAB求解
许多离散LTI都可用如下的线性常系数的差分方程描述
（1-151）
其中x[k]、y[k]分别系统的输入和输出。
在已知差分方程的N个初始状态y[k],
和输入x[k],就可由下式迭代计算出系统的输出。
y[k]=-
利用MATLAB提供的filter函数，可方便地计算出上述差分方程的零状态响应。
filter函数调用形式为
y=filter（b,a,x）
其中
a=[a0,a1,…,aN],b=[bo,b1,…,bM]
分别表示差分方程系数。
X表示输入序列，y表示输出序列。
输出序列的长度和序列相同。
例1-40试用M=9点滑动平均系统
y[k]=
处理例1-39中的受噪声干扰的正弦信号。
解：
由式（1-151）可知，M点滑动平均系统可看成是N=0的差分方程。
在调用filter函数时，调用参数a=1,b为有M个元素的向量，b中每个元素的值均为1/M。
%program1_4滑动平均去噪
M=9;
N=100;k=0:
N;
s=10*sin（0.02*pi*k）;
x=s+randn（1,N+1）;
b=ones（M,1）/M;a=1;
y=filter（b,a,x）;
plot（k,y,s,’:
’）;
xlabel（‘k’）;
ylabel（‘幅度’）
legend（‘y[k]’,’s[k]’）;
程序的运行结果如图1-52所示。
图中的虚线表示未受噪声干扰的原信号s[k],实线为9点滑动的响应y[k]。
比较图1-50的信号可以看出，系统滤出了受干扰信号中的大部分噪声，输出信号相对输入信号有4个样本的延迟。
1.9.4利用MATLAB计算DTET
当序列的DTET可写成
的有理多项式时，可用MATLAB信号处理工具箱提供的freqz函数计算DTFT的抽样值。
另外，可用MATLAB提供的abs、angle、real、imag等基本函数计算DTFT的幅度、相位、实部、虚部。
若X（
）可表示为
X（
）=
=
则freqz的调用形式为
X=freqz（b,a,w）（1-153）
式（1-153）中的b和a分别是表示式（1-152）中分子多项式和分母多项式系数的向量，即
b=[b0,b1,…,bM]
a=[a0,a1,…,aN]
w为抽样的频率点，在以式（1-153）形式调用freqz函数时，向量w的长度至少为2。
返回值X就是DTFT在抽样点w上的值。
注意一般情况下，函数freqz的返回值X是复数。
例1-41已知离散系统的H（z）为
H（z）=
试画出该系统的幅度响应。
解：
%program1_5计算系统幅度响应
b1=[0.5009-1.00190.5009];
b2=[0.53201.06400.5320];
a1=[1.0000-0.85190.4167];
a2=[1.00000.85190.4167];
b=conv（b1,b2）;%计算H（z）的分子多项式
a=conv（a1,a2）;%计算H（z）的分母多项式
w=linspace（0,pi,512）;
H=freqz（b,a,w）;
plot（w/pi,abs（H））;
ylabel（‘幅度’）;
xlabel（‘Normalizedfrequency’）;
系统幅度响应如图1-53
1.9.5部分分式法的MATLAB实现
MATLAB的信号处理工具箱提供了一个计算X（z）的部分分式展开的函数，它的调用形式如下
[r,p,k]=resifuez（b,a）
其中调用参数b,a分别表示用z
表示X（z）的分子和分母多项式。
如果X（z）的部分分式展开为
X（z）=
则residuez的返回参数r,p,k分别为
r=[r
r
r
r
]
p=[p1p2p3p3]
k=[k1k2]
同一极点p3在向量p中出现了两次，表示p3是一个二阶的重极点。
Residuez也用于由r,p,k计算z
表示的X（z）的分子和分母多项式，其调用形式为
[b,a]=residuez（r,p,k）
例1-43试用MATLAB计算
X（z）=
的部分分式展开。
解：
计算部分分式展开的MATLAB程序如下
%program1_6部分分式展开
b=[1.5,0.98,-2.608,1.2,-0.144];
a=[1,-1.4,0.6,-0.072];
[r,p,k]=residuez（b,a）;
disp（‘留书’）；disp（r’）
disp（‘极点’）；disp（p’）
disp（‘常数’）；disp（k’）;
程序运行结果为
留数
极点
常数
12
部分分式展开的结果为
X（z）=
z反变换的结果为
）u[k]+2
[k-1]
1.9.6用MATLAB计算系统的零极点
MATLAB信号处理工具箱提供的tf2zp和zp2tf函数可进行系统函数的不同表示形式的转换。
z正幂有理多项式表示的系统函数为
（1-154）
用零点、极点和常数表示的一阶因子形式的系统函数为
（1-155）
MATLAB函数[z,p,k]=tf2zp（b,a）将有理多项式表示的系统函数转换为一阶因子形式的系统函数。
[b,a]=zp2tf（z,p,k）将一阶因子形式的系统函数转换为有理多项式表示的系统函数。
例1-44试将下面的系统函数表示为一阶因子形式。
（1-156）
解：
用z正幂有理多项式表示的系统函数为
%program1_7确定一阶因子形式的系统函数
b=[100.040];
a=[1-0.80.16-0.128];
[z,p,k]=tf2zp（b,a）;
disp（‘零点’）;disp（z’）;
disp（‘极点‘）；disp（p’）;
disp（‘常数’）；disp（k’）;
程序运行结果为
零点
极点
常数
1
系统函数的一阶因子形式为
（1-157）
为了在H（z）的表达式中不出现复数，也可将式（1-157）等价的写成二阶因子的形式
（1-158）
当H（z）是表示为式（1-154）的形式时，可用[z,p,k]=tf2zp（b,a）求出系统的零极点，从而可根据极点的位置判断系统的稳定性，也可用函数zplane（b,a）画出z平面的零极点分布来判断系统的稳定性。
例１-45 已知一因果的ＬＴＩ系统的函数为
试用函数zplane画出系统的零极点分布，并判断系统的稳定性。
解：
%program1_8系统零极点分布
b=[121];
a=[1-0.5-0.0050.3];
zplane（b,a）;
图１－５４画出了系统零极点分布。
图中符号
表示零点。
符号
旁的数字表示零点的阶数，符号x表示极点，图中的虚线画的是单位圆。
由图可知，该系统的极点全在单位圆内，故系统是稳定的。
１.９.７ 简单数字滤波器的设计
例1-46 已知信号x[k]中含有角频率分别为
和
的离散正弦信号。
试设计一高通滤波器，滤出信号x[k]中低频分量，保留信号x[k]高频分量。
解：
为了简单起见，假设数字滤波器是一个具有如下形式的长度为３的ＦＩＲ系统
h[0]=h[2]=
h[1]=
系统的查分方程
+
系统的频率响应为
由上式可知系统的群延迟为
为了滤除低频分量，系统需满足
＝
可用ＭＡＴＬＡＢ命令Ａ／Ｂ求解上面的的方程组得：
下面的MATLAB程序实现了滤波器的设计及信号的滤波。
％progam1_9简单数字滤波器的设计
％确定滤波器系数
w1=0.015*pi;W2=0.4*pi;N=50;
A=[2*cos（w1）1;2*cos（w1）1];B=[0;1]
x=A/B;
a=1;b=[x
（1）x
（2）x
（1）];
%产生两个余弦信号
k=0:
N;
x1=cos（w1*k）;
x2=cos（w2*k）;
％信号滤波
y=filter（b,a,x1+x2）;
plot（k,y,`r`,k,x2,`b--`,x1+x2,`k:
`）;
ylabel（`幅度`）;xlabel（`k`）
legend（`y[k]`,`x2[k]`,`x1[k]+x2[k]`）;
图１－５５画出了程序运行结果。
由图可以看出在k=0，１时，输出响应有波动。
这是因为系统响应中存在瞬态响应。
当
时，系统进入稳态响应。
滤波后的信号相对与原信号有一个样本的延迟，着是因为在
时，
１.９.８ 语音信号的滤波
例1-47 已知一段语音信号中混入了一频率f
=1200Hz正弦型干扰信号。
语音信号的抽样频率f
=22050Hz。
试用二阶带阻滤波器滤除语音信号中的噪音信号。
解：
干扰信号的数字频率
为
由式（1-110）可得
取
，由式（1-111）可得
解上述方程得
。
当
的值为1.0619时，系统有一个极点在单位圆外，当
的值为0.941７时，系统的二个极点在单位圆内。
为了保证系统的稳定，取
＝0.9417。
由
和
的值可得系统函数为
%program1_10语音信号滤波
Fn=1200;w
=0.06;
%读入语音信号
［x,Fs］=wavread（`sample`）;
%播放语音信号
sound（x,Fs）;pause;
%设计滤波器
w0=2*pi*Fn/Fs;
beta=cos（w0）;
alpha=min（roots（[1-2/cos（w
）1]））;
a=[1,-beta*（1+alpha）,alpha];
b=[1-2*beta1]*（1+alpha）/2;
%信号滤波
y=filter（b,a,x）;
%播放滤波后语音信号
sound（y,Fs）;
图１－５６画出了滤波前后语音信号的频谱。
滤波前的语音信号的频谱在１２００Ｈz有一高峰，这是干扰所造成的。
由滤波后语音信号的频谱可看出，滤波器成功地制造了语音信号中的干扰。
2.6利用MATLAB实现信号DFT的计算
2.6.1利用MATLAB计算信号DFT
在MATLAB信号处理工具箱中，函数dftmtx（N）可用来产生N
N的DFT矩正D
。
N
N的IDFT矩正D
可用函数conj（dfmtx（N））/N来确定。
此外，MATLAB提供了4个内部函数用于计算DFT和IDFT，它们分别是：
fft（x）,fft（x,N）,ifft（X）,ifft（X,N）
fft（x）计算M点的DFT。
M是序列x的长度，即M=length（x）。
fft（x,N）计算N点的DFT。
若M>N，则将原序列截短为N点序列，再计算其N点DFT；若Mifft（X）计算M点的IDFT。
M是序列X的长度。
ifft（X,N）计算N点IDFT。
若M>N，则将原序列截短为N点序列，再计算其N点IDFT；若M为了提高fft与ifft的运算效率，应尽量使序列长度M=2
，或对序列补零使N=2
。
实现例2-7的程序如下：
%program2_1计算有限长余弦信号频谱
N=30；%数据长度
L=512；%DFT的点数
f1=100;f2=120;fs=600;
T=1/fs;
ws=2*pi*fs;
t=（0:
N-1）*T;
f=cos（2*pi*f1*t）+cos（2*pi*f2*t）;
F=fftshift（fft（f,L））;
w=（-ws/2+（0:
L-1）*ws/L）/（2*pi）;
plot（w,abs（F））;
ylabel（`幅度谱`）
实现例2-8的程序为：
%program2_2利用Hamming计算有限长余弦信号频谱
N=50；%数据长度
L=512；%DFT的点数
f1=100;f2=120;fs=600;
%N=25；
T=1/fs;
ws=2*pi*fs;
t=（0:
N-1）*T;
f=cos（2*pi*f1*t）+cos（2*pi*f2*t）;
wh=（hamming（N））`;
f=f.*wh;
F=fftshift（fft（f,L））;
w=（-ws/2+（0:
L-1）*ws/L）/（2*pi）;
plot（w,abs（F））;
ylabel（`幅度谱`）
例2-11已知一长度为16点的有限序列
试用MATLAB计算序列x[k]的16点和512点DFT。
解：
%progam2_3NumericalComputationofDTFTUsingFFT
k=0:
15;
f=cos（2*pi*k*4./16）;
F_16=fft（f）;
F_512=fft（f,512）;
L=0:
511;
plot（L/512,abs（F_512））;
holdon;
plot（k/16,abs（F_16）,`0`）;
set（gca,`xtick`,[0,0.25,0.5,0.75,1]）;
set（gca,`ytick`,[0,2,4,6,8]）;
gridon;
xlabel（`Normalizedfrequency`）;
ylabel（`Magnirude`）;
holdoff
从图2-19可见，对长度为16的序列x[k]进行512点DFT，可以获得频谱函数X（
）更多的细节，因为序列N点的离散傅立叶变换
就是序列X（
）一个周期的N个等间隔抽样，尽管从信息表达的角度，由序列x[k]16点DFTX[m]足以恢复原序列x[k]。
利用MATLAB实现由DFT计算线性卷积
例2-12试利用MATLAB实现由DFT计算有限序列线性卷积，并与线性卷积直接计算的结果进行比较。
解：
在MALAB中，序列线性卷积的直接结果是通过函数conv（x,h）来实现的。
%program2_4linearConvolutionthroughDFT
x=[1201];
h=[2211];
%Determinethelengthoftheresultofconvolution
L=length（x）+length（h）-1;
%ComputetheDFTbyzero-padding
XE=fft（x,L）;
HE=fft（h,L）;
%DeterminetheIDFToftheproduct
y1=ifft（XE.*HE）;
%plotthesequencegeneratedbycircularconvolution
%andtheerrorfromdirectlinearconvolution
k=0:
L-1;
subplot（2,1,1）;
setm（k,real（y1））;axis（[0607]）;
title（`ResultofLinearConvolution`）;
xlabel（`Timeindexk`）;ylabel（`Amplitude`）;
y2=conv（x,h）;
error=y1-y2;
subplot（2,1,2）;
stem（k,abs（error））;
xlabel（`Timeindexk`）;ylabel（`Amplitude`）;
title（`ErrorMagnitude`）;
由图2-20（b）可见，由DFT计算的线性卷积的误差非常小，这些误差主要是计算机在计算DFT时，由计算机有限字长而产生的舍入误差与计算误差。
此外，MATLAB信号处理工具箱中提供了fftfilt函数，该函数是基于重叠相加法的原理，可以实现长序列与短序列之间的线性卷积，其调用形式为
y=fftfilt（h,x,n）
其中
h:
短序列；
x：
长序列；
n：
DFT的点数。
一般选择n为2正整次幂，以提高DFT运算的效率。
3.6线性调频z变换算法
例3-3已知序列x[k]=
（
）,试计算序列x[k]在单位圆上的CZT，其中
。
解：
由于是计算序列x[k]在单位圆上的CZT，故A0=1，W0=1。
计算序列czt的MATLAB程序如下：
%program3_1
%ComputeCZTofsequencex[k]
N=13;%lengthofx[k]
n=0:
N-1;
xn=（0.8）.^n;
sita=pi/4;%phaseofstartpoint
fai=pi/6;%phaseinterval
A=exp（j*sita）;%complexstartingpoint
W=exp（-j*fai）;%complexratiobetweenpoints
M=8;%lengthofczt
y=czt（xn,M,W,A）;%callcztfunction
subplot（2,1,1）
stem（n,xn）;
xlabel（`k`）;title（`sequencex[k`]）;
m=0”M-1;
subplot（2,1,2）;
stem（m,abs（y））;
xlabel（`m`）;title（`CZTofx[k]`）;
运行结果如图3-22所示。
4.5用MATLAB实现滤波器设计
4.5.1用MATLAB实现模拟低通滤波器的设计
MATLAB信号处理工具箱提供了常用的设计模拟低通滤波器的MATLAB函数。
在实现应用中，可方便地调用这些函数完成模拟滤波器的设计。
关于这些函数现分别介绍如下：
Butterworth滤波器
[N，wc]=buttord（wp,ws,Ap,As,`s`）
[num,den]=butter（N,wc,`s`）
根据BW型滤波器的设计指标，利用MATLAB函数buttord获得BW型滤波器参数N和wc。
函数buttord的输入参数wp和ws（rad/s）分别表示滤波器的通带和阻带截频，Ap和As（dB）表示滤波器的通带和阻带衰减。
`s`表示所设计的是模拟滤波器。
函数buttord返回参数N为BW滤波器的阶数，wc（rad/s）等于BW滤波器3dB截频
。
由于wc由阻带方程确定，故由参数N、wc得出的滤波器在阻带刚好满足设计指标，在通带将超过设计指标。
当BW型滤波器的参数N、wc确定后，可用MATLAB函数butter获得BW型滤波器系统函数
的分子多项式（num）和分母多项式（den）。
ChebyshevI型滤波器
[N，wc]=cheblord（wp,ws,Ap,As,`s`）
[num,den]=cheby1（N,Ap,wc,`s`）
MATLAB函数cheblord返回参数N表示CBI型滤波器的阶数，wc（rad/s）等于wp。
参数N、wc的取值可使cheby1设计出的CBI型滤波器在通带刚好满足设计指标。
cheby1函数利用参数N、wc和Ap确定CBI型滤波器系统函数
的分子多项式（num）和分母多项式（den）。
ChebyshevII型滤波器
[N，wc]=cheb2ord（wp,ws,Ap,As,`s`）
[num,den]=cheby2（N,As,wc,`s`）
MATLAB函数cheb2ord返回参数N表示CBII型滤波器的阶数，wc（rad/s）的取值可使cheby2设计出的滤波器在通带刚好满足设计指标。
cheby2函数利用参数N、wc和As确定CBII滤波器系统函数
的分子多项式（num）和分母多项式（den）。
椭圆滤波器
[N，wc]=ellipord（wp,ws,Ap,As,`s`）
[num,den]=ellip（N,Ap,wc,`s`）
MATLAB函数ellipord返回参数N表示椭圆滤波器的阶数，wc=wp。
MATLAB函数ellip确定阶数为N，通带衰减为Ap（dB）,阻带衰减为As（dB）,通带截频为wc的椭圆滤波器系统函数
的分子多项式和分母多项式。
在滤
升级会员 