DSP上机matlab函数使用详细解答56页Word文件下载.docx
《DSP上机matlab函数使用详细解答56页Word文件下载.docx》由会员分享,可在线阅读,更多相关《DSP上机matlab函数使用详细解答56页Word文件下载.docx(43页珍藏版)》请在冰点文库上搜索。
![DSP上机matlab函数使用详细解答56页Word文件下载.docx](https://file1.bingdoc.com/fileroot1/2023-5/9/6a2bd509-52ba-4a3c-ba95-ab2f93273530/6a2bd509-52ba-4a3c-ba95-ab2f932735301.gif)
其中调用参数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的取值范围
若不加此句则当做kx和kh从0开始
stem(k,y);
xlabel(‘k’);
ylabel(‘y’);
程序的运行结果如图1-51所示。
1.9.3离散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对应y,b对应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;
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,’:
’);
%改成(k,y,k,s,‘:
’用虚线表示)不能用stem
xlabel(‘k’);
ylabel(‘幅度’)
legend(‘y[k]’,’s[k]’);
%图标说明
程序的运行结果如图1-52所示。
图中的虚线表示未受噪声干扰的原信号s[k],实线为9点滑动的响应y[k]。
比较图1-50的信号可以看出,系统滤出了受干扰信号中的大部分噪声,输出信号相对输入信号有4个样本的延迟。
1.9.4利用MATLAB计算DTFT
当序列的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);
%卷积计算系数!
!
a=conv(a1,a2);
%计算H(z)的分母多项式
w=linspace(0,pi,512);
%默认是512点,常常在0到pi间
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
]
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’);
程序运行结果为
留数
0.70000.50000.3000
极点
0.60000.60000.2000
常数
12
部分分式展开的结果为
z反变换的结果为
x[k]=(0.7x0.6
+0.5x(k+1)0.6
+0.3x0.2
)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];
%此处的0千万别少,少了阶数就少了,一定要考虑到z正幂的常数项即z负幂的最低次项以及所缺的项!
a=[1-0.80.16-0.128];
[z,p,k]=tf2zp(b,a);
disp(‘零点’);
disp(z’);
%此处有无分号都会显示
disp(‘极点‘);
disp(p’);
disp(‘常数’);
零点
00+0.2000i0-0.2000i
0.80000.0000+0.4000i0.0000-0.4000i
1
系统函数的一阶因子形式为
(1-157)
为了在H(z)的表达式中不出现复数,也可将式(1-157)等价的写成二阶因子的形式
(1-158)
当H(z)是表示为式(1-154)的形式时,可用[z,p,k]=tf2zp(b,a)求出系统的零极点,从而可根据极点的位置判断系统的稳定性,也可用函数zplane(b,a)(用于给出参数直接根据154式子画图!
)画出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]=
系统的查分方程
+
系统的频率响应为
由上式可知系统的群延迟为看exp(-j
)
为了滤除低频分量,系统需满足
=
保留omga2
可用MATLAB命令A/B求解上面的的方程组得:
下面的MATLAB程序实现了滤波器的设计及信号的滤波。
%progam1_9简单数字滤波器的设计
%确定滤波器系数
w1=0.015*pi;
w2=0.4*pi;
N=50;
A=[2*cos(w1)1;
2*cos(w2)1];
B=[0;
1]%分号表示换行
x=A\B;
%B左乘一个A的逆即A\
a=1;
b=[x
(1)x
(2)x
(1)];
%取X中元素用X(),方向为以列优先取
%产生两个余弦信号
x1=cos(w1*k);
要先构造关于x的载波才能使用filter函数
x2=cos(w2*k);
%信号滤波
y=filter(b,a,x1+x2);
%为何为x1+x2
plot(k,y,`r`,k,x2,`b--`,k,x1+x2,`k:
`);
%冒号表示用虚线
ylabel(`幅度`);
xlabel(`k`)
legend(`y[k]`,`x2[k]`,`x1[k]+x2[k]`);
%分别对应1,2,3
图1-55画出了程序运行结果。
由图可以看出在k=0,1时,输出响应有波动。
这是因为系统响应中存在瞬态响应。
当
时,系统进入稳态响应。
滤波后的信号相对与原信号有一个样本的延迟,着是因为在
时,
1.9.8 语音信号的滤波(带通滤波器怎样设计的,下面那几个公式出处)
例1-47 已知一段语音信号中混入了一频率f
=1200Hz正弦型干扰信号。
语音信号的抽样频率f
=22050Hz。
试用二阶带阻滤波器滤除语音信号中的噪音信号。
干扰信号的数字频率
为
%即W数字角频率
由式(1-110)可得
取
,由式(1-111)可得
解上述方程得
的值为1.0619和0.9417。
的值为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)先对x离散才能DFT变换。
fft(x,N)计算N点的DFT。
若M>
N,则将原序列截短为N点序列,再计算其N点DFT;
若M<
N,则将原序列补零至N点,然后计算其N点DFT。
ifft(X)计算M点的IDFT。
M是序列X的长度。
ifft(X,N)计算N点IDFT。
N,则将原序列截短为N点序列,再计算其N点IDFT;
N,则将原序列补零至N点,然后计算其N点IDFT。
为了提高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;
%对cos函数以Ts抽样
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);
%决定ws的映射方式,也可以直接写ws*(0:
L-1)
plot(w,abs(F));
ylabel(`幅度谱`)
实现例2-8的程序为:
%program2_2利用Hamming计算有限长余弦信号频谱
N=50;
L=512;
f1=100;
%N=25;
wh=(hamming(N))`;
f=f.*wh;
%用汉明实现就在是使f乘以以汉明系数
F=fftshift(fft(f,L));
例2-11已知一长度为16点的有限序列
试用MATLAB计算序列x[k]的16点和512点DFT。
N=16;
r=4;
L=512;
ws=2*pi*r/N;
k=0:
N-1;
f=cos(ws*k);
F=fft(f,L);
w=ws*(0:
ylabel('
幅度谱'
)我自己的。
%progam2_3NumericalComputationofDTFTUsingFFT
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));
L长度和F长度要一样才行
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;
ylabel(‘Magnirude’);
holdoff
从图2-19可见,对长度为16的序列x[k]进行512点DFT,可以获得频谱函数X(
)更多的细节,因为序列N点的离散傅立叶变换
就是序列X(
)一个周期的N个等间隔抽样,尽管从信息表达的角度,由序列x[k]16点DFTX[m]足以恢复原序列x[k]。
2.6.2利用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
L-1;
%?
?
subplotreal是啥子
subplot(2,1,1);
%(subplot(m,n,p)将输出窗口分为m*n矩阵,p是当前窗口的标号,标号大小从按行优先数)
stem(k,real(y1));
axis([0607]);
%说横坐标为0到6纵为0到7,若果不规定在幅度到6的时候会截断
title(`ResultofLinearConvolution`);
xlabel(`Timeindexk`);
ylabel(`Amplitude`);
y2=conv(x,h);
%直接计算线性卷积,x,h为序列,切都从0开始
error=y1-y2;
subplot(2,1,2);
stem(k,abs(error));
%为啥子要加个abs喃error也是时域的
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:
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%调用czt计算z变换,xn许留函数,之后为点数,fai和sita的那个负数形式
subplot(2,1,1)
stem(n,xn);
xlabel(`k`);
title(`sequencex[k]`);
m=0:
M-1;
stem(m,abs(y));
xlabel(`m`);
title(`CZTofx[k]`);
运行结果如图3-22所示。
4.5用MATLAB实现滤波器设计
4.5.1用MATLAB实现模拟低通滤波器的设计
MATLAB信号处理工具箱提供了常用的设计模拟低通滤波器的MATLAB函数。
在实现应用中,可方便地调用这些函数完成模拟滤波器的设计。
关于这些函数现分别介绍如下:
Butterwort