MATLAB程序设计语言在信号处理中的应用.docx
《MATLAB程序设计语言在信号处理中的应用.docx》由会员分享,可在线阅读,更多相关《MATLAB程序设计语言在信号处理中的应用.docx(32页珍藏版)》请在冰豆网上搜索。
MATLAB程序设计语言在信号处理中的应用
1概述
1.1MATLAB程序设计语言简介
MATLAB,MatrixLaboratory的缩写,是由Mathworks公司开发的一套用于科学工程计算的可视化高性能语言,具有强大的矩阵运算能力。
与大家常用的Fortran和C等高级语言相比,MATLAB的语法规则更简单,更贴近人的思维方式,被称之为“草稿纸式的语言”。
截至目前,MATLAB已经发展到12.1版,适用于所有32位的Windows操作系统,按NTFS(NT文件系统)格式下完全安装约需850MB。
MATLAB软件主要由主包、仿真系统和工具箱三大部分组成。
1.2MATLAB应用入门
1.MATLAB的安装与卸载
MATLAB软件在用户接口设计上具有较强的亲和力,其安装过程比较典型,直接运行光盘中的安装向导支撑程序SETUP.exe,按其提示一步步选择即可。
MATLAB自身带有卸载程序,在其安装目录下有uninstall子目录,运行该目录下的uninstall.exe即可;也可以通过Windows系统的安装卸载程序进行卸载。
2.MATLAB的启动与退出
MATLAB安装完成后,会自动在Windows桌面上生成一个快捷方式,它是指向安装目录下\bin\win32\matlab.exe的链接,双击它即可来到MATLAB集成环境的基本窗口,通常称之为命令窗口。
MATLAB的退出与普通WIN32的程序一样,值得一提的是它有一个自身专有的快捷键Ctrl+Q。
3.MATLAB界面简介
图1MATLAB基本界面——命令窗口
1)菜单栏
菜单栏中包括File、Edit、View、Web、Window和Help六个菜单项。
这里着重介绍File项。
File项是数据输入/输出的接口,包括10个子项,这里重点介绍其中的5个子项:
New:
新建文件项。
有四个选择:
MFile(.M,文本格式的MATLAB程序文件,可以直接通过文件名的方式在MATLAB环境下解释运行);Figure(图形);Model(仿真模型文件)和GUI(可视化界面文件)。
Open:
打开所有MATLAB支持的文件格式,系统将自动识别并采用相应的程序对文件进行处理。
例如,打开一个.m文件,系统将自动打开M文件编辑器对它进行编辑。
ImportData...:
导入用于MATLAB处理的数据函数,包括各种图像文件、声音文件和.mat文件。
SaveWorkspaceAs...:
将工作空间的变量以.mat(二进制)或ASCII文本的形式存入文件。
SetPath...:
设置工作路径。
可以打开路径设置(SetPath)对话框(图8-2),将用户自己建立的目录加入MATLAB的目录系统中,以便所编制的文件能够在MATLAB环境中直接调用。
图2路径设置对话框
单击AddFolder...按钮可以将你的一个文件夹加入到系统路径中;AddwithSubfolders...允许把一个文件夹包括其所有的子文件夹加入到系统路径中。
这两种操作均可以直观地在右侧的路径栏内看到结果。
选中一个加入的文件夹,你可以利用MovetoTop(移至所有路径的最前面),MoveUp(上移一个),MoveDown(下移一个),MovetoBottom(移至所有路径的最后面)等四个按钮将改变文件在系统路径中的排列位置以利于对文件的搜索使用,也可以利用Remove按钮将其删除。
对路径操作完毕后,按Save按钮予以保存;按Close按钮关闭本对话框;按Revert按钮取消所有未保存的改动;按Default按钮将还原到MATLAB安装时的路径设置;按Help按钮则启动帮助系统解答疑难。
2)命令行区
对输入命令的解释MATLAB按以下顺序进行:
①检查它是否是工作空间中的变量,是则显示变量内容。
②检查它是否是嵌入函数,是则运行之。
③检查它是否是子函数。
④检查它是否是私有函数。
⑤检查它是否是位于MATLAB搜索路径范围内的函数文件或脚本文件。
请注意,如果有两个以上的方案与输入的命令相匹配,MATLAB将只执行第一个匹配。
4.MATLAB常用命令
2基本数值运算
2.1MATLAB内部特殊变量和常数
MATLAB内部有很多变量和常数,用以表达特殊含义。
常用的有:
(1)变量ans:
指示当前未定义变量名的答案。
(2)常数eps:
表示浮点相对精度,其值是从1.0到下一个最大浮点数之间的差值。
该变量值作为一些MATLAB函数计算的相对浮点精度,按IEEE标准,eps=2-52, 近似为2.2204e-016。
(3)常数Inf:
表示无穷大。
当输入或计算中有除以0时产生Inf。
(4)虚数单位i,j:
表示复数虚部单位,相当于。
(5)NaN:
表示不定型值,是由0/0运算产生的。
(6)常数pi:
表示圆周率π,其值为3.1415926535897…。
2.2变量类型
1.变量命名规则
MATLAB中对变量的命名应遵循以下规则:
(1)变量名可以由字母、 数字和下划线混合组成,但必须以字母开头。
(2)字符长度不能大于31。
(3)变量命名区分大小写。
2.局部变量和全局变量
局部变量是指那些每个函数体内自己定义的,不能从其他函数和MATLAB工作空间访问的变量。
全局变量是指用关键字“global”声明的变量。
全局变量名应尽量大写,并能反映它本身的含义。
如果需要在工作空间和几个函数中都能访问一个全局变量,必须在工作空间和这几个函数中都声明该变量是全局的。
2.3矩阵及其运算
MATLAB具有强大的矩阵运算和数据处理功能,对矩阵的处理必须遵从代数规则。
1.矩阵生成
1)一般矩阵的生成
对于一般的矩阵MATLAB的生成方法有多种。
最简单的方法是从键盘直接输入矩阵元素。
直接输入矩阵元素时应注意:
各元素之间用空格或逗号隔开,用分号或回车结束矩阵行,用中括号把矩阵所有元素括起来。
例1在工作空间产生一个3×3矩阵A可用MATLAB语言描述如下:
A=[123;456;789]
或A=[123
456
789]
运行结果:
A=
123
456
789
2)特殊矩阵的生成
对于特殊的矩阵可直接调用MATLAB的函数生成。
用函数zeros生成全0矩阵:
格式B=zeros(m,n)生成m×n的全0阵。
用函数ones生成全1矩阵:
格式B=ones(m,n)生成m×n的全1阵。
用函数eye生成单位阵:
格式B=eye(m,n)生成m×n矩阵,其中对角线元素全为1,其他元素为0。
2.矩阵的运算
矩阵的运算有基本运算和函数运算两种类型。
基本运算包括矩阵的加、减、乘、除、乘方、求转置、求逆等,其主要特点是通过MATLAB提供的基本运算符+、-、*、/(\)、^等即可完成。
函数运算主要是通过调用MATLAB系统内置的运算函数来求取矩阵的行列式(det(A)),求秩(rank(A)),求特征值和特征向量([V,D]=eig(A)),求Jordan标准形(jordan(A))和矩阵分解等。
需要用时可以参阅联机帮助和相关参考书。
例2矩阵的基本运算。
A=[1,2,3;4,5,6];
B=[6,5,4;3,2,1];
C=A+B%计算两个矩阵的和
D=B′%计算矩阵B的转置
E=A*D%做矩阵乘法,必须要满足矩阵乘法的基本要求
%E应该是2阶方阵
F=det(E)%求E的行列式值
G=E^(-1)%求E的逆
输出结果:
C=
777
777
D=
63
52
41
E=
2810
7328
F=54
G=
0.5185-0.1852
-1.35190.5185
3基本语句
3.1程序控制语句
1. 循环语句
MATLAB的循环语句包括for循环和while循环两种类型。
1)for循环
语法格式:
for循环变量=起始值:
步长:
终止值
循环体
end
起始值和终止值为一整形数,步长可以为整数或小数,省略步长时,默认步长为1。
执行for循环时,判定循环变量的值是否大于(步长为负时则判定是否小于)终止值,不大于(步长为负时则小于)则执行循环体,执行完毕后加上步长,大于(步长为负时则小于)终止值后退出循环。
例3给矩阵A、B赋值。
MATLAB语句及运行结果如下:
k=5;
a=zeros(k,k)%矩阵赋零初值
form=1∶k
forn=1∶k
a(m,n)=1/(m+n-1);
end
end
fori=m∶-1∶1
b(i)=i;
end
运行结果:
a=
1.00000.50000.33330.25000.2000
0.50000.33330.25000.20000.1667
0.33330.25000.20000.16670.1429
0.25000.20000.16670.14290.1250
0.20000.16670.14290.12500.1111
b=
12345
2)while循环
语法格式:
while表达式
循环体
end
其执行方式为:
若表达式为真(运算值非0),则执行循环体;若表达式为假(运算结果为0),则退出循环体,执行end后的语句。
例4a=3;
whilea
a=a-1
end
输出:
a=2
a=1
a=0
2.条件转移语句
条件转移语句有if和switch两种。
1)if语句
MATLAB中if语句的用法与其他高级语言相类似,其基本语法格式有以下几种:
格式一:
if逻辑表达式
执行语句
end
格式二:
if逻辑表达式
执行语句1
else
执行语句2
end
格式三:
if逻辑表达式1
执行语句1
elseif逻辑表达式2
执行语句2
end
2)switch语句
switch语句的用法与其他高级语言相类似,其基本语法格式为:
switch表达式(标量或字符串)
case值1
语句1
case值2
语句2
…
otherwise
语句n
end
3.2绘图语句
常用的MATLAB绘图语句有figure、plot、subplot、stem等,图形修饰语句有title、axis、text等。
1.figure
figure有两种用法,只用一句figure命令,会创建一个新的图形窗口,并返回一个整数型的窗口编号。
figure(n)表示将第n号图形窗口作为当前的图形窗口,并将其显示在所有窗口的最前面;如果该图形窗口不存在,则新建一个窗口,并赋以编号n。
2.plot
线型绘图函数。
用法为plot(x,y,′s′)。
参数x为横轴变量,y为纵轴变量,s用以控制图形的基本特征如颜色、粗细等,通常可以省略,常用方法如表2所示。
3.Stem
绘制离散序列图,常用格式stem(y)和stem(x,y)分别和相应的plot函数的绘图规则相同,只是用stem命令绘制的是离散序列图。
4.subplot
subplot(m,n,i)图形显示时分割窗口命令,把一个图形窗口分为m行,n列,m×n个小窗口,并指定第i个小窗口为当前窗口。
5.绘图修饰命令
在绘制图形时,我们通常需要为图形添加各种注记以增加可读性。
在plot语句后使用title(′标题′)可以在图形上方添加标题,使用xlabel(′标记′)或ylabel(′标记′)为X轴或Y轴添加说明,使用text(X值、Y值、′想加的标示′)可以在图形中任意位置添加标示。
例5画图基本语句如图3所示。
MATLAB语句及运行结果如下:
x=0:
0.1*pi:
2*pi;%定义x向量
figure
(1);%创建一个新的图形窗口,编号为1
subplot(2,2,1);%将窗口划分为2行,2列,在第1个窗口中作图
plot(x,sin(x));%画图
title(′正弦线′);%给图形加标题
subplot(2,2,2);%在第2个窗口中作图
plot(x,sin(x),′r′);%画一正弦波,红色
xlabel(′X′);%给x轴加说明
ylabel(′SIN(X)′);%给y轴加说明
subplot(2,2,3);%在第2个窗口中作图
plot(x,sin(x),′--′);%画一正弦波,破折线
subplot(2,2,4);%在第2个窗口中作图
plot(x,sin(x),′r+′);%画一正弦波,红色破折线
text(4,0,′注记′);
图3例5中绘制的几种正弦波形
4MATLAB函数
4.1函数及其调用方法
在MATLAB语言中,M文件有两种形式:
脚本和函数。
脚本没有输入/输出参数,只是一些函数和命令的组合。
它可以在MATLAB环境下直接执行,也可以访问存在于整个工作空间内的数据。
由脚本建立的变量在脚本执行完后仍将保留在工作空间中可以继续对其进行操作,直到使用clear命令对其清除为止。
函数是MATLAB语言的重要组成部分。
MATLAB提供的各种工具箱中的M文件几乎都是以函数的形式给出的。
函数接收输入参数,返回输出参数,且只能访问该函数本身工作空间中的变量,从命令窗或其他函数中不能对其工作空间的变量进行访问。
1.函数结构
MATLAB语言中提供的函数通常由以下五个部分组成:
(1)函数定义行;
(2)H1行;
(3)函数帮助文件;
(4)函数体;
(5)注释。
这五个部分中最重要的是函数定义行和函数体。
函数定义行:
MATLAB语言在M文件的第一行用关键字“function”把M文件定义为一个函数,并指定它的名字(必须和文件名相同),同时也定义了函数的输入和输出参数。
函数定义行是一个MATLAB函数所必需的,其他各部分的内容可以没有,这种函数称为空函数。
例如:
求最大值函数“max”的定义行可描述为
function[Y,I]=max(x)
其中,“max”为函数名,输入参数为“x”,输出参数为“Y”和“I”。
函数体:
函数体是函数的主体部分,它包括进行运算和赋值的所有MATLAB程序代码。
函数体中可以包括流程控制、输入/输出、计算、赋值、注释以及函数调用和脚本文件调用等。
在函数体中完成对输出参数的计算。
2.函数调用
函数调用的过程实际上就是参数传递的过程。
例如,在一个脚本文件里调用函数“max”可采用如下方式:
n=1:
20;
a=sin(2*pi*n/20);
[Y,I]=max(a);
该调用过程把变量“a”传给了函数中的输入参数“x”,然后把函数运算的返回值传给输出参数“Y”和“I”。
其中,Y是a序列的最大值,I是最大值Y对应的坐标值。
4.2常用数字信号处理函数
1.信号产生函数
1)三角波或锯齿波发生函数:
sawtooth()
语法格式:
sawtooth(t,width)。
产生以2π为周期幅值范围在[-1,+1]之间的三角波或锯齿波。
参数t为时间向量;width是[0,1]之间的数,它决定函数在一个周期内上升部分和下降部分的比例。
width=0.5产生三角波,width=1产生锯齿波,此时函数可简写为:
sawtooth(t)。
2)方波发生函数:
square()
语法格式:
square(t)。
产生以2π为周期幅值范围在[-1,+1]之间的方波,参数t为时间向量。
3)sinc发生函数:
sinc()
语法格式:
sinc(t)
例6信号产生举例
clearall
t=0:
0.0001:
0.1;
x1=sawtooth(2*pi*50*t);%在[0,0.1]之间产生5个周期的锯齿波
subplot(221)
plot(t,x1)
x2=sawtooth(2*pi*50*t,0.5);%在[0,0.1]之间产生5个周期的三角波
subplot(222)
plot(t,x2)
x3=square(2*pi*50*t);%在[0,2]之间产生10个周期的方波
subplot(223)
plot(t,x3)
axis([0,0.1,-1.2,1.2])
t=-4:
0.1:
4;
x4=sinc(t);%产生抽样函数
subplot(224)
plot(t,x4)
运行结果如图4所示。
图4常用信号
(a)锯齿波;(b)三角波;(c)方波;(d)抽样函数
2.常用窗的MATLAB函数表示
表3常用窗的MATLAB函数表示
窗名称
MATLAB函数
窗名称
MATLAB函数
矩形窗
boxcar(N)
哈明窗
hamming(N)
三角窗
triang(N)
布莱克曼窗
blackman(N)
汉宁窗
hanning(N)
凯塞-贝尔窗
kaiser(N,BETA)
说明:
除凯塞-贝尔窗外其他窗函数的使用方法相同。
函数的参数N是窗长度,调用结果为一个列向量。
例产生50点的哈明窗可用MATLAB语言表示为:
y=hamming(50);plot(y)凯塞-贝尔窗函数是一组可调窗函数。
其语法格式为:
Kaiser(N,BETA),返回一个N点的Kaiser窗,参数BETA是窗函数表达式中的参数β,其含义参照前面的理论部分介绍。
3.滤波器分析与实现函数
1)取绝对值:
abs()
语法格式:
abs(x)。
当x为实数时计算x的绝对值;x为复数时得到的是复数的模值;x为字符串时得到各字符的ASCII码。
2)取相角:
angle()
语法格式:
angle(z)。
求复矢量或复矩阵的相角,结果为一个以弧度为单位介于-π和+π之间的值。
3)求线性卷积:
conv()
语法格式:
conv(x,y)。
求矢量x和y的卷积,若x(n)和y(n)的长度分别为M和N,则返回值是长度为M+N-1的矢量。
例7x(n)=[345];y(n)=[2678],求其线性卷积。
MATLAB语句如下:
x=[345];
y=[2678];
z=conv(x,y)
运行结果:
z=
62655826740
4)利用指定的数字滤波器对数据进行滤波:
filter()
常用语法格式:
y=filter(b,a,x)。
函数filter利用数字滤波器对数据进行滤波时,采用直接Ⅱ型结构实现,因而适用于IIR和FIR两种滤波器。
参数:
a=[a0a1a2…aM],b=[b0b1b2…bN]是滤波器系数,x为输入序列矢量,y为滤波后的输出。
即:
滤波器的系统函数为:
标准形式中取a0=1,若输入滤波系数a中a0≠1时,MATLAB会自动归一化系数;若a0=0,系统给出出错信息。
例8在语音信号处理中,常利用周期脉冲信号通过AR(10)模型来近似合成浊音信号。
若信号周期T=46,AR模型的系数a=[1,-1.7218,1.2594,-0.6157,0.7754,-0.6496,0.3651,-0.4547,0.3339,0.0975,-0.1851],试合成5个周期的浊音信号。
用MATLAB语句实现如下:
T=46;
a=[1,-1.7218,1.2594,-0.6157,0.7754,-0.6496,0.3651,-0.4547, 0.3339, 0.0975,-0.1851];fori=0:
5
x(i*T+1)=1;
end
y=filter(1,a,x);
plot(y)
图5合成浊音波形
5)计算数字滤波器H(z)的频率响应H(ejω):
freqz()
语法格式:
[H,W]=freqz(B,A,N)得到数字滤波器的N点的频率向量W和与之相对应的N点的频率响应向量H,计算所得的N个频率点均匀的分布在[0,π]上。
参数A=[a0a1a2…aM],B=[b0b1b2…bN]是滤波器系数,即滤波器H(z)形式如下:
参数N(与上式中阶次N的含义不同)最好选用2的整数次幂,以便使用FFT进行快速运算,N的缺省值为512。
freqz(B,A,N)直接绘制频率响应图,而不返回任何值。
H=freqz(B,A,W)返回W向量中指定的频率范围内的频率响应。
其中,W以弧度为单位在[0,π]范围内。
[H,F]=freqz(B,A,N,Fs)对H(ejω)在[0,Fs/2]上等间隔采样N点,采样点频率及相应的频率响应值分别记录在F和H中。
6)计算数字滤波器H(z)的单位脉冲响应h(n):
impz()
语法格式:
[H,T]=impz(B,A)。
滤波器用传递函数模型限定,参数B,A分别为H(z)分子分母多项式的系数,函数返回滤波器的冲击响应列向量H和时间即采样间隔列向量T。
4.变换函数
1)一维快速离散Fourier变换:
fft()
语法格式:
y=fft(x)。
y是计算信号x的快速离散傅里叶变换。
当x为矩阵时,计算x中每一列信号的离散傅里叶变换。
当x的长度为2的幂时,用基2算法;否则,采用较慢的分裂基算法。
y=fft(x,n)。
计算n点的FFT。
当x的长度大于n时,截断x;否则补零。
2)一维快速离散Fourier逆变换:
ifft()
语法格式:
y=ifft(x)。
y是计算信号x的快速离散傅里叶变换的逆变换。
y=ifft(x,n)。
计算n点的快速离散傅里叶变换的逆变换。
3)离散余弦变换(DCT):
dct()
语法格式:
y=dct(x)。
计算信号x的离散余弦变换。
y=dct(x,n)。
计算n点的离散余弦变换。
当x的长度大于n时,截断x;否则补零。
离散余弦逆变换可由函数idct实现。
例9计算信号
n=1,2,…,100的DCT。
用MATLAB语言可实现如下:
N=100;n=1:
N;
x=n+20*sin(2*pi*n/20);
y=dct(x);
z=id