MATLAB在时间序列分析中的应用

《MATLAB语言》课程论文

MATLAB在时间序列分析中的应用

姓名：李伟

学号：12010245258

专业：通信工程

班级：10级通信一班

指导老师：汤全武

学院：物理电气信息学院

完成日期：2011.12.10

MATLAB 在时间序列分析中的应用

（李伟 12010245258 2010级通信1班）

[摘要]Matlab 强大的科学计算和可视化功能使其在各个领域中得到了广泛的应用．采用Maflab 进行时序列分析可以极大地简化编程工作，并具有界面友好、操作方便的特点．介绍了使用Matlab 进行时间序列分析的基本方法和步骤，并通过实例进行了说明．

[关键词]Matlab 时间序列 谱分析

一、问题的提出

1984年美国的MathWorks 公司推出了Mat —lab ，在许多领域得到了充分的利用．其强大的科学计算与可视化功能，开放式的可扩展环境以及其各种功能强大的工具箱(ToolBox)，使得它成为计算机辅助设计与分析、算法研究和应用开发的基本工具和首选平台．时间序列分析是采用参数模型对观测得到的有序随机数据进行分析的一种处理方法，通过时间序列可以对系统的动态特性进行分析、对系统的状态进行预测，从而为系统的状态监控和故障诊断提供依据．Matlab 工具箱中包含了许多函数，借助于这些函数可以方便地实现系统的时间序列分析．

二、时间序列分析原理及实现

时间序列分析(autoRegressive movingAverage)是对有序的随机数据(信号)处理的一种方法，它的出发点是承认数据的有序性和相关性，通过数据内部的相互关系来辨识系统的变化规律，它的建模方法是将系统的输出看作是在白噪声输入下的响应．具体地讲，就是针对一组试验数据，建立系统的参数模型，ARMA(m ，n)的参数模型可以表示为：

10()()()m n

x t x t t τλτλ?τ?σλ==+-=-∑∑ (1) 式中：{ (t)， (t 一1)， (t 一2)? 。 (t —m)}为有序的时间序列，{ (t)， (t 一1)，? ， (t —A)}为有序的白噪声序列，方程的左边为系统的自回归部分，它反映了系统的固有特性，右边表示系统的滑动平均部分，当0τ?=时为MA 模型，当0λ?=时为AR 模型．辨识系统模型参数的方法有很多种，常用的方法主要有最dx-乘法、辅助变量法、Marple 法等．根据不同的需要和研究对象可以采用不同的建模方法．在建立了系统的模型后，可以对系统的状态进行预测、分析预测误差、进行谱分析．关于这些算法的基本原理，可以参考文献[2～4]，这些在Matlab 中都提供了相应的函数．

采用Matlab 进行时间序列分析主要包括4步．

1)数据的读人

Matlab 采用类似于C 语言的方式进行数据的读人，可以直接从数据文件中将数据读到一个矩阵中．

fid=fopen(fileName ，，r ，)； ％打开一个文件进行读写

data=fscanf(fid ，'g')； ％ 将数据读人到data 中

status=fclose(fid)； ％ 释放文件句柄

2)建立模型

在获得所要分析的数据后可以对数据进行建模，本文主要介绍2个函数： th = ar(y ，n)；

h = ivar(y ，n)；

ar(y ，n)函数采用最小二乘法进行模型参数的估计，该函数要求输入噪声为白噪声，当输入噪声为色噪声时，不能保证模型参数的估计值的无偏性和一致收敛性．而ivar(y ，n)则采用最优辅助变量的方法进行参数的估计，计算得到的参数模型存放在th 中．th 中的数据采用Matlab 独有的THETA 格式模型进行定义．通过th2arx()函数可以得到模型参数和THETA 格式的转换．

3)模型分析

模型的分析包括模型的仿真、预测及误差分析和谱分析．

e= pe(y ，th)；

y1= idsim(y ，th)；

y1= predict(y ，th)；

y2=th2ff(th)；

pe(y ，th)用于计算模型实测值与估计值之间的误差，误差值存放在e 中. idsim (y ，th)对输入的数据进行仿真，并将仿真结果存放在y 中.dict(y ，th)则针对模型的输人数据和模型格式进行预测，并将预测值存放在Y 中，th2f(th)可以实现求数据的频响函数．

4)图形输出

Matlab 提供了强大的数据输入输出的功能．对数据的分析结果，可以采用图形的方式进行直观的表示，常用的针对时间序列分析的绘图函数有： Plot(x,y1,x,y2 )，在同一个图中对分析结果进行表示．

Bodeplot(e)，直接画出波德图．

Ffplot(e)，画频谱图．

Nyqplot()，画奈氏图．

三、Theta 模型参数

Theta 格式是Matlab 系统辨识工具箱中通用的参数模型格式，Theta 模型的定义可以分为两种，即基于输入输出表示的Theta 模型和基于状态空间表示的Theta 模型．基于输入输出的Theta 模型可以对应各种输入输出参数模型，如AR 、ARX 、ARMA 、BJ 等；基于状态空间表示的Theta 模型则与连续或离散状态空间参数模型对应，它们的信息都以矩阵的形式存储，但模型信息数据的组织结构不同．在时间序列分析中，常用的是第一种数据模型，其结构可以表示为：

1111()()()()()()()()()()()

n nu nu n B q B q C q A q y t u t nk u t nk e t F q F q D q =-+Λ+-+ （2） 公式(2)中，A(q)、B(q)、F(q)( =1，2? ．^n)、C(q)、D(q)为平移算子q 的多

项式，其阶次分别为M 、nbi、nfi( =1，2?．，n)、ItC、nd，nix为系统的输入变量个数．设n为所有多项式的阶次之和，令r=max(凡，7，6+3nu)，则系统的输入输出The—ta模型格式为如下定义的(3+凡)X r矩阵，矩阵中每行的内容表示为：

1)矩阵的第1行为估计方差，采样时间，输入个数眦，各个多项式的阶次M 、nb、ItC、nd、nf,nk；

2)第2行为最终预测误差FPE，模型生成的13期、时间和命令；

3)第3行为估计参数的向量，即A、、C、D、F 的系数；

4)第4行到第3+凡行为估计的方差矩阵；

5)对于连续系统，该矩阵可能增加到凡+4行，其中包含系统的死区时间．对于时间序列分析而言，在生成Theta模型以后，需要根据不同的需要对该模型进行分析，以便从中提取所需的估计参数以及最终的误差．

四、应用实例

为了对上面的方法进行说明，采用5个问题加以说明，第1个问题利用小波时间序列进行消噪或压缩.第2个问题利用ddencmp和wdencmp函数实现数据降噪.第3个问题利用函数wavedec对时间序列进行一维多分辨分析.第4个问题用wthcoef对时间序列的小波系数进行阈值处理.第5个问题利用wprcoef，由wpdec 得到的t对时间蓄力分解的一位小波包系数重构.5个问题的源程序如下：

问题1

格式:

（1）[xc,cxc,lxc,perf0,petfl2]

=wdencmp('gbl,x,'wname',N,thr,sorh,keepapp)

（2）[xc,cxc,lxc,perf0,perfl2]=wdencmp('lvd',x,'wname',N,thr,sorh) （3）[xc,cxc,lxc,perf0,pefgl2]=wdencmp('lvd',C,L,'wname',N,thr,sorh) 说明:

x：待消噪或压缩的时间序列。

N：小波分解层数。

c,l：由函数wavedec得到的时间序列一维小波变换的分解系数及其维数。

'wname'：消噪时选取的小波函数。

thr：阈值。

sorh：软、硬阈值的选择。

keepapp：可选参数，当keepapp=1时表示低频系数不进行量化，反之低频系数要进行量化。

xc：发表会消噪或压缩后的时间序列。

cxc,lxc：xc的小波分解结构。

perf0,perfl2：恢复和压缩L平方的范数百分比。

'gbl'：global的缩写，各层用同一个阈值。

'lvd'：level-dependent的缩写，没一层用不同的阈值。

[xc,cxc,lxc,perf0,perfl2]=('gbl',X,'wname

',N,thr,sorh,keepapp) %全局正阈值消噪或压缩。

[xc,cxc,lxc,perf0,perfl2]=wdencmp('lvd',x, 'wname',N,thr,SORH)

%利用ddencmp函数返回的thr,sorh消噪或压缩。

[xc,cxc,lxc,perf0,perfl2]=wdencmp('lvd',c,l, 'wname',N,thr,sorh)

%利用wavedec函数返回的c,l消噪或压缩。

例：

load leleccum;

indx=2600:3100;

x=leleccum(indx); %原始信号、

Subplot(221),plot(x);axis([0,500,100,500]);

title('原始信号');

thr=35;

[xd,cxd,lxd,perf0,perfl2]=wdencmp('gbl',x,'db3',3,thr,'h',1);%固定阈值压缩处理

subplot(222),plot(xd);

axis([0,500,100,500]);

title('固定阈值压缩信号');

[thr,sorh,keepapp]=ddencmp('den','wv',x); %缺省值消噪或压缩

xd=wdencmp('gbl',x,'db3',2,thr,sorh,keepapp); %利用ddencmp的返回值进行消噪或压缩

subplot(223),plot(xd);

axis([0,500,100,500]);

title('利用gbl的压缩信号');

xd=wdencmp('lvd',x,'db3',2,thr,sorh); %利用ddencmp函数返回的thr,sorh消噪或压缩

Subplot(224),plot(xd);

axis([0,500,100,500]);

title('利用lvd的压缩信号');

运行结果如图1

图1 原始信号及其不同方式下的消噪或压缩信号

问题2

%由nosimima函数生成一个一维信号，使用小波分析对信号进行降噪处理

load noismima;

indx=1:1000;

x=noismima(indx);

[thr,sorh,keepapp]=ddencmp('den','wv',x); %得到降噪参数

xd=wdencmp('gbl',x,'db3',2,thr,sorh,keepapp);%用全局阈值进行降噪处理subplot(211) %画图

plot(x)

subplot(212)

plot(xd)

运行结果如图2

图2 原始序列及其降噪后的序列

问题3

格式：

（1）[c,1]=wacedec(x,N,’wname ’)

(2) [c,1]=wacedec(x,N,Lo_D,Hi_D)

说明：

x ：待分析的时间序列。

N ：分解层数。

’wname ’：选取的小波函数。

Lo_D,Hi_D ：分别为小波分解滤波器的低频系数和高频系数。

c ：时间序列的以为多分辨分解系数11[,,,......,]j j j c cA cD cD cD -=，其中j cA 第j 阶分解的低频系数，11,,......,j j cD cD cD -分别为第j,j-1,……,1阶的高频系数。

L ：11,,,......,j j j cA cD cD cD -对应的序列长度组成的序列。

例：

load sumsin;

s=sumsin; %调用数据sumsin

[c,l]=wavedec(s,3,'db1'); %利用db1分解3阶

subplot(211);

plot(s);

title('原始信号');\

subplot(212);

plot(c);

title('小波分解系数');

运行结果如图3

图3 原始信号解气小波分解系数

问题4

格式：

（1）NC=wthcoef('d',c,l,N,P)

（2）NC=wthcoef('d',c,l,N)

（3）NC=wthcoef('a',c,l)

（4）NC=wthcoef('t',c'l'N,T,SORH)

说明：

c,l:由函数wavedec得到的时间序列以为小波变换的分解系数及其维数。 N:包含高频尺度的向量。

P：把较小小波系数置零的百分比信息向量。

sorh：软，硬阈值的选择。

T：小波树结构。

NC：阈值处理后的[NC,l]即为[c,l]经过阈值处理后的结果。

NC=wthcoef('d',c,l,N,P)%对小波分解结构[c,l]进行压缩，把较小小波新书全部置零。

NC=wthcoef('d',c,l)%对小波分解结构[c,l]进行压缩，把高频小波系数全部置零，返回分解结构NC。

NC=wthcoef('a',c,l)%对小波分解结构[c,l]进行压缩，把低频小波系数全部置零，返回分解结构NC。

NC=wthcoef('t',c,l,N,t,sorh)%对小波分解结构[c,l]进行阈值压缩处理，当sorh='s'时进行软阈值处理；当sorh='h'时进行硬阈值处理。

例：

load leleccum;

indx=2600:3100;

x=leleccum(indx);

Subplot(311),plot(x);axis([0,500,100,500]);

title('原始信号');

[c,l]=wavedec(x,3,'db3');

N=[1,2,3];%设置尺度向量

P=[98,99,97];%设置阈值向量

NC=wthcoef('d',c,l,N,P,'s');%对高频小波系数进行软阈处理，把高频小波系数全部置零

x1=waverec(NC,l,'db3'); %重构

subplot(312),plot(x1);axis([0,500,100,500]);title('软阈值重构信号'); NC=wthcoef('t',c,l,N,P,'h');%对高频小波系数进行硬阈值处理，把高频小波系数全部置零

x2=waverec(NC,l,'db3'); %重构

subplot(313),plot(x1);axis([0,500,100,500]);title('硬阈值重构信号'); 运行结果如图4：

图4 原始信号及其不同方式下阈值处理后的重构信号

问题5

格式：

（1)x=wprcoef(t,N)

(2)x=wprcoef(t)

说明：

t：分解后的小波树结构。

N：小波树的节点。

x：时间序列的一维小波包系数重构。

例：

load noisdopp;

x=noisdopp; %调用信号

subplot(121);

plot(x);title('原始信号');

t=wpdec(x,3,'db3','shannon');

rcfs=wprcoef(t,[2,1]); %重构（2,1）上信号

Subplot(122);

plot(rcfs);

title('小波包节点（2,1）处重构信号');

运行结果如图5

图5 原始信号及其节点（2,1）上的重构信号

五、结论

通过上面的实例可以看出，利用Madab工具，可以很方便地实现数据的时间序列建模，建模精度高、分析结果可靠，同时利用它强大的图形界面功能可以很直观地进行数据的分析与比较．当然，并不是所有时间序列的问题都可以用Madab 内带的函数实现，有些需要用户自己编制响应的程序．但是，借助于它的强大的计算功能，可以很方便地实现这些功能．

六、课程体会

通过一学期对Matlab语言的学习，我收获甚多.在相较别的计算机语言，我

深刻体会到Matlab语言简单易懂，对我们各门课程有非常大的帮助.Matlab是一个高级矩阵/阵列预言，它包含控制语句、函数、数据结构、输入和输出和面向对象编程的特点。用可以在命令窗口中将输入语句与执行命令同步，也可以先编好一个较大的复杂的应用程(M文件)后再一起欲行。新版本的MATLAB预言师基于最为流行的C++预言基础上的，因此语法特征与C++预言极为相似，而且更加简单，更加符合科技人员对数学表达式的书写格式。使之更利于计算机专业的科技人员使用。而且这种语言的可移植性好、课拓展性极强，这也是MATLAB能够深入到科学研究及工程计算各个领域的重要原因。

在科学研究和工程应用中，为了客服一般预言对大量数学运算，尤其当涉及矩阵运算时编制程序麻烦等困难，应运而生了MATLAB编程运算的软件，在自动控制、图像处理、语言处理、信号分析、振动理论、优化设计、时序分析和系统建模等领域都能得到很好的处理效果。而且在MATLAB中，可以直接在Simulink环境中运作的工具包很多，已覆盖通信、控制、信号出题、DSP、电力系统等诸多领域，所涉及的内容专业性极强。

参考文献：

[1] 徐听．Matlab工具箱应用指南——控制工程篇[M]．北京：人民邮电出版社，2000．

[2] 黄世霖．工程信号处理[M]．北京：人民交通出版社，1986．

[3] 温熙森．机械系统动态分析理论与应用[M]．长沙：国防科技大学出版社，1998．

[4] 黄仁．机械设备工况监测与故障诊断[M]．南京：东南大学出版社，1990．