6. 小波分析 wavelet analysis

自从学习过佛瑞艾尔变形和频率估计，作者对小波分析产生的兴趣，开始阅读一些相关的资料。同时发现正在给自己上数学课的David Walnut 教授是小波分析的泰山北斗。他写的这方面的教课书遍布全球。作者的心中又产生的无比的崇敬与羡慕，所以将当前世界最先进的小波分析技巧写出来与大家分享。

读过《R语言时间系列中文教程》都应该知道如何使用弗瑞艾尔变形估计频率。但必须假设，被估计的频率是始终存在于波动之中的。更经常的状况是在一整个波中某一频率只在这个波中的一小部分出现。使用弗瑞艾尔变形不可能监测到在某一时间点上的频率变化，因为它假设所估计的频率都是自始至终存在的。例如，下面的波中有一个很慢的频率是始终存在的，在中间部分突然出现了频率非常高的新波动，而且很快就消失了。这样的波动需要使用小波分析。

t=1:500

c1 = 2*cos(2*pi*t/150 + .6*pi)

plot.ts(c1)

t2= ifelse(t>200 & t<300,t,0)

c2=1*cos(2*pi*t2/10 + .6*pi)

par(new=T)

plot.ts(c2)

par(new=F)

例如心电图是有来测量人心脏跳动的手段，心脏的各个组成部分不是同时不停的工作的，在一个心房收缩的过程中某些的心房是休息的。所以在心电图上来看，微小的波动是突然出现仅仅持续很短的时间就消失了，像这样的微小波动就要使用小波分析来捕捉。

心电图是医生进行临床诊断的方法之一，通过阅读心电图医生可以推测病人的心脏是正常的还是哪里出现的毛病。接下来我们要介绍使用小波分析计算机智能进行自动诊断，也就是说我们积累了很多人的心电图数据通过小波分析将心电图的特色提取出来。这些特色将告诉我们病人的心脏是正常的还是有状况。通过这些心电图的数据我们可以建立一个数学模型，当有一个新的病人进来我们就可以对他的心电图进行诊断。

这里使用的数据也是网上的数据，有600行，每一行代表着一个波动，虽然这些波动不是心电图的波动，但我们把它全当做心电图的波动来使用，只是为了介绍概念。(A) Downward Trend. (B) Cyclic. (C) Normal. (D) Upward Shift. (E) Upward Trend. (F) Downward Shift.

从上面的贴图可以看到波动分为6种，第一种是向下型，第二种是循环型，第三种是正常型，第四种是向上移动型，第五种是向上型，第六种是向下移动型。

在数据中，1到100 行为第一型，101 到200 行为第二型，以此类推。

# extracting DWT coefficients (with Haar filter)

library(wavelets)

feature<-NULL

mydata <- read.table("http://archive.ics.uci.edu/ml/databases/synthetic_control/synthetic_control.data",header=F, sep="")

#mydata <- read.table("C://Users//User//Desktop//R Language//Wavelet//synthetic_control.data",header=F, sep="")

for (i in 1:nrow(mydata)) {

a <- t(mydata[i,])

wt <- dwt(a, filter="haar", boundary="periodic")

feature <- rbind(feature, unlist(c(wt@W,wt@V[[wt@level]])))

}

feature <- as.data.frame(feature)

上面的命令是用来读取数据提取小波分析数据特色的，所使用的程序包叫做WAVELETS。 数据是通过使用read.table 命令直接读取一个因特网的连接。上面的FOR 循环是从1循环到600 ，也就是对于每一行的数据都要执行FOR 循环的命令。其中 最关键的命令为DWT 命令（离散小波），这个命令把每一行波动进行小波分析并且提出其中的特色。其中我们使用的是其中HAAR 小波，BOUNDARY 的设置为PERIODIC 也就是循环的。RBIND 命令只是把所有的特色困绑在一起存储于FEATURE 变量中。最后一句的命令是将生成的FEATURE矩阵转化为DATA.FRAME 对象可做下面的使用。这里是上面一部分代码运行的贴图

# set class labels into categorical values

classId <- c(rep("1",100), rep("2",100), rep("3",100),

rep("4",100), rep("5",100), rep("6",100))

wtSc <- data.frame(cbind(classId, feature))

我们数据中有6种波，上面的命令是将波的类型（1-6）和波的特色捆绑在一起生成新的DATA.FRAME对象称作wtSc 。

# build a decision tree with ctree() in package party

library(party)

ct <- ctree(classId ~ ., data=wtSc,

controls = ctree_control(minsplit=30, minbucket=10, maxdepth=5))

pClassId <- predict(ct)

# check predicted classes against original class labels

table(classId, pClassId)

上面一节的语句中使用了一个叫作PARTY 的程序包，作用是在于对前面提取出的FEATURE 对象建立决策树模型。CTREE 命令是用来生成这个决策树的，来看一下这一节输出结果

输出中有一个6 X 6 的矩阵被称作“对错矩阵”。可以看到我们原本有100个1型的波动，其中有97 个被这个数学模型正确的认识出来了，有3 个被错误的认作2型波动。再比方说，这个数学模型正确的认出了99个2型波动，但其中的一个被错误的认作1型波动。

# accuracy

(sum(classId==pClassId)) / nrow(wtSc)

这条语句用来计算生成的决策树模型的正确率，值为87%。即为百分之八十七的正确率。

plot(ct, ip_args=list(pval=FALSE), ep_args=list(digits=0))

PLOT 语句用来输出决策树模型的插图

比方说，数据第一行代表的病人，或者是一个与这个病人状况非常相似的病人，回到这里从新做了心电图。数据中1-100行都是1型波动。通过我们的数学模型，还可以知道这个病人的心电图是1型的吗？通过数学模型来判定病人的心电图型号被称作预测。下面就是进行预测的代码，同样需要将心电图的特色提取出来，然后放回到数学模型中来判定型号。

#make predictions

f<-NULL 初始化F 变量

a <- t(mydata[1,]) 1号病人的心电图存入a 中，代表这个人又回来了。

wt <- dwt(a, filter="haar", boundary="periodic") 输入提取特色

f <- rbind(f, unlist(c(wt@W,wt@V[[wt@level]]))) 特色放入f 中

f <- as.data.frame(f) f 转为DATA.FRAME对象

predict(ct,f) 使用前面模型的名字叫ct, 特色叫f 来进行预测

预测结果是 [1] 1 是说明我们的数学模型已经认出这个病人的心电图是1 型的。