MATLAB程序设计
杨凯 2010 . 11
内容回顾 一、二维绘图
关键命令关键命令::plot 二、三维绘图
关键命令关键命令:: 三维曲线三维曲线::plot3 三维曲面
三维曲面::meshgird;mesh;surf
主要内容——数据分析与计算 一、数据统计处理 二、数据插值 三、曲线拟合 四、离散傅立叶变换 五、多项式计算
一、数据统计处理 1.1最大值和最小值
MATLAB 提供的求数据序列的最大值和最小值的函数分别为max 和min ,两个函数的调用格式和操作过程类似数的调用格式和操作过程类似。。
1、求向量的最大值和最小值
求一个向量X 的最大值的函数有两种调用格式格式,,分别是分别是::
(1 y=max(X:返回向量X 的最大值存入y ,如果X 中包含复数元素中包含复数元素,,则按模取最大值则按模取最大值。。
(2 [y,I]=max(X:返回向量X 的最大值存入y ,最大值的序号存入I ,如果X 中包含复数元素中包含复数元素,,则按模取最大值取最大值。。
求向量X 的最小值的函数是min(X,用法和max(X完全相同完全相同。。 例:求向量x 的最大值的最大值。。 命令如下命令如下:: x=[-43,72,9,16,23,47];
y=max(x %求向量x 中的最大值
[y,l]=max(x %求向量x 中的最大值及其该元素的位置
2.求矩阵的最大值和最小值
求矩阵A 的最大值的函数有3种调用格式种调用格式,,分别是分别是::
(1 max(A:返回一个行向量返回一个行向量,,向量的第i 个元素是矩阵A 的第i 列上的最大值列上的最大值。。
(2 [Y,U]=max(A:返回行向量Y 和U ,Y 向量记录A 的每列的最大值的每列的最大值,,U 向量记录每列最大值的行号大值的行号。。
(3 max(A,[],dim:dim取1或2。dim取1时,该函数 ;dim取2时,该函数返回一个和max(A完全相同 完全相同; 行上的最大值。 。
,其第i个元素是A矩阵的第i行上的最大值列向量, 列向量 完全相同。 。
求最小值的函数是min,其用法和max完全相同 例:分别求3×4矩阵x中各列和各行元素中的最大 。
值,并求整个矩阵的最大值和最小值
并求整个矩阵的最大值和最小值。
3.两个向量或矩阵对应元素的比较
函数max 和min 还能对两个同型的向量或矩阵进行比较进行比较,,调用格式为调用格式为::
(1 U=max(A,B:A,B 是两个同型的向量或矩阵是两个同型的向量或矩阵,,结果U 是与A,B 同型的向量或矩阵同型的向量或矩阵,,U 的每个元素等于A,B 对应元素的较大者对应元素的较大者。。
(2 U=max(A,n:n 是一个标量是一个标量,,结果U 是与A 同型的向量或矩阵的向量或矩阵,,U 的每个元素等于A 对应元素和n 中的较大者中的较大者。。
min 函数的用法和max 完全相同完全相同。。
例:求两个2×3矩阵x, y 所有同一位置上的较大元素构成的新矩阵p 。
1.2 求和与求积
数据序列求和与求积的函数是sum和
是一个向量,,其使用方法类似。。设X是一个向量prod,其使用方法类似 是一个矩阵,,函数的调用格式为
函数的调用格式为:: A是一个矩阵 各元素的和。。
sum(X:返回向量X各元素的和 各元素的乘积。。
prod(X:返回向量X各元素的乘积
返回一个行向量,,其第i个元素是A sum(A:返回一个行向量 。
的第i列的元素和 列的元素和。
prod(A:返回一个行向量返回一个行向量,,其第i 个元素是A 的第i 列的元素乘积的元素乘积。。
sum(A,dim:当dim 为1时,该函数等同于sum(A;当dim 为2时,返回一个列向量返回一个列向量,,其第i 个元素是A 的第i 行的各元素之和行的各元素之和。。
prod(A,dim:当dim 为1时,该函数等同于prod(A;当dim 为2时,返回一个列向量返回一个列向量,,其第i 个元素是A 的第i 行的各元素乘积行的各元素乘积。。
例: 求矩阵A 的每行元素的乘积和全部元素的乘积的每行元素的乘积和全部元素的乘积。。
1.3 平均值和中值
求数据序列平均值的函数是mean ,求数据序列中值的函数是median 。两个函数的调用格式为两个函数的调用格式为::mean(X:返回向量X 的算术平均值的算术平均值。。
median(X:返回向量X 的中值的中值。。
mean(A:返回一个行向量返回一个行向量,,其第i 个元素是A 的第I 列的算术平均值列的算术平均值。。
median(A:返回一个行向量返回一个行向量,,其第i 个元素是A 的第i 列的中值列的中值。。
mean(A,dim:当dim 为1时,该函数等同于mean(A;当dim 为2时,返回一个列向量返回一个列向量,,其第i 个元素是A 的第I 行的算术平均值行的算术平均值。。
当为时,等同于;当median(A,dim:dim 1median(Adim 为2时,返回一个列向量返回一个列向量,,其第i 个元素是A 的第i 行的中值的中值。。
例:分别求向量x 与y 的平均值和中值的平均值和中值。。
1.4 标准方差与相关系数 1.求标准方差
在MATLAB 中,提供了计算数据序列的标准方差的函数std 。对于向量X ,std(X返回一个标准方差方差。。对于矩阵A ,std(A返回一个行向量返回一个行向量,,它的各个元素便是矩阵A 各列或各行的标准方差各列或各行的标准方差。。std 函数的一般调用格式为的一般调用格式为::
Y=std(A,flag,dim
其中dim 取1或2。当dim=1时,求各列元素的标准方差;当dim=2时,则求各行元素的标准方差则求各行元素的标准方差。。flag 取0或1
例: 对二维矩阵x ,从不同维方向求出其标准方差从不同维方向求出其标准方差。。
2、排序
MATLAB 中对向量X 是排序函数是sort(X,函数返回一个对X 中的元素按升序排列的新向量中的元素按升序排列的新向量。。
sort 函数也可以对矩阵A 的各列或各行重新排序,其调用格式为其调用格式为:: [Y,I]=sort(A,dim
其中dim 指明对A 的列还是行进行排序的列还是行进行排序。。若dim=1,则按列排则按列排;;若dim=2,则按行排则按行排。。Y 是排序后的矩阵后的矩阵,,而I 记录Y 中的元素在A 中位置中位置。。例:对二维矩阵做各种排序对二维矩阵做各种排序。。
自学
其他部分统计函数((自学 3.其他部分统计函数 (1累计和cumsum (2累积积sumprod (3相关系数corrcoef
这些函数调用的参数与操作方式都与MEDIAN 因此不作详细的介绍。 。
,因此不作详细的介绍函数基本上一样 (中值 中值
函数基本上一样,
二、数据插值
插值的定义——是对某些集合给定的数据点之间函数的估值方法数的估值方法。。
当不能很快地求出所需中间点的函数时当不能很快地求出所需中间点的函数时,,插值是一个非常有价值的工具个非常有价值的工具。。
2.1 一维数据插值
在MATLAB 中,实现这些插值的函数是interp1,其调用格式为其调用格式为:: Y1=interp1(X,Y,X1,'method'
函数根据X,Y 的值的值,,计算函数在X1处的值处的值。。X,Y 是两个等长的已知向量是两个等长的已知向量,,分别描述采样点和样本值分别描述采样点和样本值。。
Y1=interp1(X,Y,X1,'method'
X1是一个向量或标量是一个向量或标量,,描述欲插值的点描述欲插值的点,,Y1是一个与X1等长的插值结果等长的插值结果。。method 是插值方法是插值方法,,允许的取值有许的取值有‘‘linear’、‘nearest’、‘cubic’、‘spline’。
注意注意::X1的取值范围不能超出X 的给定范围的给定范围,,否则,会给出会给出““NaN”错误错误。。
MATLAB 中有一个专门的3次样条插值函数
Y1=spline(X,Y,X1,其功能及使用方法与函数Y1=interp1(X,Y,X1,‘spline’完全相同完全相同。。
例:插值在车轮外形曲线检测中的应用车轮曲线特点车轮曲线特点::多段圆弧状曲线拼接而成
基本原理基本原理::磁爬式
系统组成系统组成::(1两个旋转编码器两个旋转编码器;; (2磁性小球磁性小球;;
θ1
θ2O 1O 2O 3(3机械支架机械支架;;
测量曲线:(测量曲线
例:某观测站测得某日6:00时至18:00时之间每隔2小时的室内外温度(℃,用3次样条插值分别求得该日室内外6:30至17:30时之间每隔2小时各点的近似温度(℃。
h =6:2:18;
t=[18,20,22,25,30,28,24;15,19,24,28,34,32,30]';设时间变量h 为一行向量为一行向量,,温度变量t 为一个两列矩阵矩阵,,其中第一列存放室内温度其中第一列存放室内温度,,第二列储存室外温度温度。。命令如下命令如下::
XI =6.5:2:17.5
YI=interp1(h,t,XI,‘spline’ %用3次样条插值计算
2.2 二维数据插值
在MATLAB 中,提供了解决二维插值问题的函数interp2,其调用格式为其调用格式为::
Z1=interp2(X,Y,Z,X1,Y1,'method'
其中X,Y 是两个向量是两个向量,,分别描述两个参数的采样点样点,,Z 是与参数采样点对应的函数值是与参数采样点对应的函数值,,X1,Y1是两个向量或标量个向量或标量,,描述欲插值的点描述欲插值的点。。Z1是根据相应的插值方法得到的插值结果插值方法得到的插值结果。。method 的取值与一维插值函数相同值函数相同。。X,Y,Z 也可以是矩阵形式也可以是矩阵形式。。
同样同样,,X1,Y1的取值范围不能超出X,Y 的给定范围范围,,否则否则,,会给出会给出““NaN”错误错误。。
三、曲线拟合
在MATLAB 中,用polyfit 函数来求得最小二乘拟合多项式的系数拟合多项式的系数,,再用polyval 函数按所得的多项式计算所给出的点上的函数近似值式计算所给出的点上的函数近似值。。
polyfit 函数的调用格式为函数的调用格式为:: [P,S]=polyfit(X,Y,m
函数根据采样点X 和采样点函数值Y ,产生一个m 次多项式P 及其在采样点的误差向量S 。其中X,Y 是两个等长的向量是两个等长的向量,,P 是一个长度为m+1的向量的向量,,P 的元素为多项式系数的元素为多项式系数。。
polyval 函数的功能是按多项式的系数计算x 点多项式的值
polyval 函数用来求代数多项式的值数多项式的值,,其调用格式为: Y=polyval(P,x
若x 为一数值为一数值,,则求多项式在该点的值则求多项式在该点的值;;若x 为向量或矩阵为向量或矩阵,,则对向量或矩阵中的每个元素求其多项式的值其多项式的值。。
例:已知数据表[t,y],试求2次拟合多项式p(t,然后求ti=1,1.5,2,2.5,…,9.5,10各点的函数近似值各点的函数近似值。。
例: x0=0:0.1:1;
y0=[-.447 1.978 3.11 5.25 5.02 4.66 4.01 4.58 3.45 5.35 9.22]; p=polyfit(x0,y0,3
p = 56.6915 -87.1174 40.0070 -0.9043 xx=0:0.01:1;yy=polyval(p,xx; plot(xx,yy,'-b',x0,y0,'or'
电力机车: 非自给性牵引方式 ~
输电线 牵引变电所 分区亭 馈电线接触网 电力机车 钢轨
检测原理:力传递(杠杆
已在多个机务段投入使用
回顾:如果对车轮外形曲线采用拟合? 拟合和插值使用场合:
(1拟合适用于可以用曲线方程来表示;
(2插值适用于无法用曲线方程来表示或为较复杂的分段函数的分段函数。。
四、离散傅立叶变换
自学::数字信号处理 4.1离散傅立叶变换算法简要(自学 的实现((快速傅立叶变换 4.2 离散傅立叶变换的实现 其调用格式与功能为: :
,其调用格式与功能为 函数,
一维离散傅立叶变换函数
(1 fft(X:返回向量X的离散傅立叶变换。设X的长度(即元素个数为N,若N为2的幂次,则为以2为基数的快速傅立叶变换,否则为运算速度很慢的非2幂次
。
应用于矩阵的每一列。
。对于矩阵X,fft(X应用于矩阵的每一列的算法。 的算法
(2 fft(X,N:计算N 点离散傅立叶变换。它限定向量的长度为N ,若X 的长度小于N ,则不足部分补上零;若大于N ,则删去超出N 的那些元素些元素。。对于矩阵X ,它同样应用于矩阵的每一列它同样应用于矩阵的每一列,,只是限定了向量的长度为N 。
(3 fft(X,[],dim或fft(X,N,dim:这是对于矩阵而言的函数调用格式的函数调用格式,,前者的功能与FFT(X基本相同基本相同,,而后者则与FFT(X,N基本相同基本相同。。只是当参数dim=1时,该函数作用于X 的每一列的每一列;;当dim=2时,则作用于X 的每一行的每一行。。
值得一提的是值得一提的是,,当已知给出的样本数N0不是2的幂次时次时,,可以取一个N 使它大于N0且是2的幂次,然后利用函数格式fft(X,N或fft(X,N,dim便可进行快速傅立叶变换。这样这样,,计算速度将大大加快。
相应地,相应地,一维离散傅立叶逆变换函数是ifft 。ifft(F返回F 的一维离散傅立叶逆变换;ifft(F,N为N 点逆变换;ifft(F,[],dim或ifft(F,N,dim则由N 或dim 确定逆变换的点数或操作方向的点数或操作方向。。
例:给定数学函数
x(t=12sin(2π×10t+π/4+5cos(2π×40t取N=128,试对t 从0~1秒采样秒采样,,用fft 作快速傅立叶变换速傅立叶变换,,绘制相应的振幅-频率图频率图。。在0~1秒时间范围内采样128点,从而可以确定采样周期和采样频率确定采样周期和采样频率。。由于离散傅立叶变换时的下标应是从0到N-1,故在实际应用时下标应该前移1。又考虑到对离散傅立叶变换来说换来说,,其振幅| F(k|是关于N/2对称的对称的,,故只须使k 从0到N/2即可即可。。
程序如下程序如下:: N=128; % 采样点数 T=1; % 采样时间终点
t=linspace(0,T,N; % 给出N 个采样时间
ti(I=1:Nx=12*sin(2*pi*10*t+pi/4+5*cos(2*pi*40*t; % 求各采样点样本值x
dt=t(2-t(1; % 采样周期 f=1/dt; % 采样频率(Hz
X=fft(x; % 计算x 的快速傅立叶变换X F=X(1:N/2+1; % F(k=X(k(k=1:N/2+1f=f*(0:N/2/N; % 使频率轴f 从零开始
plot(f,abs(F,'-*' % 绘制振幅-频率图 xlabel('Frequency'; ylabel('|F(k|'
例:车轮踏面擦伤动态检测系统检测系统:(:(加速度方式度方式
背景资料: (1擦伤定义 ,车轮滑 成的,
车轮擦伤是由车轮发生滑行造成的行时车轮踏面与钢轨接触的那部分成了固定的磨擦面,它与钢轨持续摩擦而使车轮踏面上发生局部平面磨耗,形成擦伤。
(2擦伤的危害
踏面擦伤超限会锤击钢轨、降低车辆运行品质,以及使车辆弹簧等零部件折损,轴承保持架裂纹或破碎,引发热轴、切轴,甚至造成列车事故。
(3检测原理
用振动加速度传感器捕捉车轮与钢轨之间的振动信号。
车轮踏面擦伤动态检测系统检测系统:(:(加速度方式度方式
课外资料课外资料::
(1《基于Matlab 中FFT 函数的电力谐波分析方法函数的电力谐波分析方法》》,
李达义,电测与仪表,2002年第7期P23-26
主要内容:对电力系统谐波电压(或电流进行了仔细分析,推导了该计算公式与FFT 函数计算出的谐波系数的关系。
课外资料课外资料::用FFT 谐波分析后的效果
课外资料课外资料::
(2《基于MATLAB 的振动信号平滑处理方法》,孙苗钟,电子测量技术,2007年6月P55—57主要内容:对振动信号进行多种平滑处理,对平滑前后的效果进行了FFT 分析。
(3《基于MATLAB 系统的信号FFT 频谱分析与显示显示》》,曾尚璀,科技通报,2000年7月P241—246主要内容:给出一种用MATLAB 系统实现信号频谱分析与显示的方法的方法,,能快速地在普通通用PC 机上完成信号的频谱分析和频谱图显示
五、多项式计算
鉴于MATLAB 无零下标下标,,故把多项式的一般形式表达为: 多项式求根 1 1
21+−++++n n n n a x a x a x a 1. 命令格式命令格式::x=roots(A。这里A 为多项式的系数A(1,A(2, ,A(N,A(N+1;解得的根赋值
给数组X ,即X(1,X(2, ,X(N。
【例*】试用ROOTS 函数求多项式x 4+8x 3-10的根
这是一个4次多项式次多项式,,它的五个系数依次为:1,8,0,0,-10。下面先产生多项式系数的向量A ,然后求根后求根::x =
A=[1800-10] A =
1800-10 x=roots(A
-8.0194 -0.5075 + 0.9736i -0.5075 -0.9736i 1.0344
2.多项式的建立
若已知多项式的全部根若已知多项式的全部根,,则可以用POLY 函数建立起该多项式该多项式;;也可以用POLY 函数求矩阵的特征多项式多项式。。POLY 函数是一个MATLAB 程序,调用它的命令格式是调用它的命令格式是::
A=poly(x 若x 为具有N 个元素的向量个元素的向量,,则poly(x建立以x 为其根的多项式为其根的多项式,,且将该多项式的系数赋值给向量A 。在此种情况下,POLY 与ROOTS 互为逆函数函数;;若x 为N ×N 的矩阵x ,则poly(x返回一个向量赋值给A ,该向量的元素为矩阵x 的特征多项式之系数多项式之系数::A(1,A(2, ,A(N,A(N+1。
