离散信号和系统分析的MATLAB实现.docx
《离散信号和系统分析的MATLAB实现.docx》由会员分享,可在线阅读,更多相关《离散信号和系统分析的MATLAB实现.docx(50页珍藏版)》请在冰豆网上搜索。
![离散信号和系统分析的MATLAB实现.docx](https://file1.bdocx.com/fileroot1/2023-2/5/11a26be0-dcde-4968-978d-538b36746454/11a26be0-dcde-4968-978d-538b367464541.gif)
离散信号和系统分析的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平面的零极点分布来判断系统的稳定性。
例1-45 已知一因果的LTI系统的函数为
试用函数zplane画出系统的零极点分布,并判断系统的稳定性。
解:
%program1_8系统零极点分布
b=[121];
a=[1-0.5-0.0050.3];
zplane(b,a);
图1-54画出了系统零极点分布。
图中符号
表示零点。
符号
旁的数字表示零点的阶数,符号x表示极点,图中的虚线画的是单位圆。
由图可知,该系统的极点全在单位圆内,故系统是稳定的。
1.9.7 简单数字滤波器的设计
例1-46 已知信号x[k]中含有角频率分别为
和
的离散正弦信号。
试设计一高通滤波器,滤出信号x[k]中低频分量,保留信号x[k]高频分量。
解:
为了简单起见,假设数字滤波器是一个具有如下形式的长度为3的FIR系统
h[0]=h[2]=
h[1]=
系统的查分方程
+
系统的频率响应为
由上式可知系统的群延迟为
为了滤除低频分量,系统需满足
=
可用MATLAB命令A/B求解上面的的方程组得:
下面的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]`);
图1-55画出了程序运行结果。
由图可以看出在k=0,1时,输出响应有波动。
这是因为系统响应中存在瞬态响应。
当
时,系统进入稳态响应。
滤波后的信号相对与原信号有一个样本的延迟,着是因为在
时,
1.9.8 语音信号的滤波
例1-47 已知一段语音信号中混入了一频率f
=1200Hz正弦型干扰信号。
语音信号的抽样频率f
=22050Hz。
试用二阶带阻滤波器滤除语音信号中的噪音信号。
解:
干扰信号的数字频率
为
由式(1-110)可得
取
,由式(1-111)可得
解上述方程得
。
当
的值为1.0619时,系统有一个极点在单位圆外,当
的值为0.9417时,系统的二个极点在单位圆内。
为了保证系统的稳定,取
=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);
图1-56画出了滤波前后语音信号的频谱。
滤波前的语音信号的频谱在1200Hz有一高峰,这是干扰所造成的。
由滤波后语音信号的频谱可看出,滤波器成功地制造了语音信号中的干扰。
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的椭圆滤波器系统函数
的分子多项式和分母多项式。
在滤