R語言小波分析wavelet

1、6.小波分析wavelet analysis自從學(xué)習(xí)過佛瑞艾爾變形和頻率估計，作者對小波分析產(chǎn)生的興趣，開始閱讀一些相關(guān)的 資料。同時發(fā)現(xiàn)正在給自己上數(shù)學(xué)課的David Walnut教授是小波分析的泰山北斗。他寫的這方面的教課書遍布全球。作者的心中乂產(chǎn)生的無比的崇敬與羨慕，所以將當(dāng)前世界最 先進(jìn)的小波分析技巧寫出來與大家分享。讀過R語言時間系列中文教程都應(yīng)該知道如何使用弗瑞艾爾變形估計頻率。但必須假 設(shè)，被估計的頻率是始終存在丁波動之中的。更經(jīng)常的狀況是在一整個波中某一頻率只在 這個波中的一小部分出現(xiàn)。使用弗瑞艾爾變形不可能監(jiān)測到在某一時間點上的頻率變化， 因為它假設(shè)所估計的頻率都是自始至終存

2、在的。例如，下面的波中有一個很慢的頻率是始 終存在的，在中間部分突然出現(xiàn)了頻率非常高的新波動，而且很快就消失了。這樣的波動 需要使用小波分析。t=1:500cl = 2*cos(2*pi*t/150 + .6*pi)plot.ts(c1)t2= ifelse(t200 & t300,t,0)c2=1*cos(2*pi*t2/10 + .6*pi)par(new=T)plot.ts(c2)par(new=F)例如心電圖是有來測量人心臟跳動的手段，心臟的各個組成部分不是同時不停的工作的， 在一個心房收縮的過程中某些的心房是休息的。所以在心電圖上來看，微小的波動是突然 出現(xiàn)僅僅持續(xù)很短的時間

3、就消失了，像這樣的微小波動就要使用小波分析來捕捉。心電圖是醫(yī)生進(jìn)行臨床診斷的方法之一，通過閱讀心電圖醫(yī)生可以推測病人的心臟是正常 的還是哪里出現(xiàn)的毛病。接下來我們要介紹使用小波分析計算機(jī)智能進(jìn)行自動診斷，也就 是說我們積累了很多人的心電圖數(shù)據(jù)通過小波分析將心電圖的特色提取出來。這些特色將 告訴我們病人的心臟是正常的還是有狀況。通過這些心電圖的數(shù)據(jù)我們可以建立一個數(shù)學(xué) 模型，當(dāng)有一個新的病人進(jìn)來我們就可以對他的心電圖進(jìn)行診斷。這里使用的數(shù)據(jù)也是網(wǎng)上的數(shù)據(jù)，有600行，每一行代表著一個波動，雖然這些波動不是 心電圖的波動，但我們把它全當(dāng)做心電圖的波動來使用，只是為了介紹概念。(A) Downwar

4、d Trend. (B) Cyclic. (C) Normal. (D) Upward Shift. (E) Upward Trend. (F)Downward Shift.從上面的貼圖可以看到波動分為6種，第一種是向下型，第二種是循環(huán)型，第 三種是正常型，第四種是向上移動型，第五種是向上型，第六種是向下移動型。在數(shù)據(jù)中，1到100行為第一型，101到200行為第二型，以此類推。# extracting DWT coefficients (with Haar filter) library(wavelets) feature-NULLmydata - read.table(http:/arch

5、/ml/databases/synthetic control/syntheticcon trol.data,header=F, sep=)#mydata - read.table(C:/Users/User/Desktop/RLanguage/Wavelet/synthetic_control.data,header=F, sep=)for (i in 1:nrow(mydata) a - t(mydatai,)wt - dwt(a, filter=haar, boundary=periodic)feature - rbind(feature, unlist(c

6、(wtW,wtVwtlevel)feature - as.data.frame(feature)上面的命令是用來讀取數(shù)據(jù)提取小波分析數(shù)據(jù)特色的，所使用的程序包叫做WAVELETS數(shù)據(jù)是通過使用read.table命令直接讀取一個因特網(wǎng)的連接。 上面 的FOR循環(huán)是從1循環(huán)到600，也就是對于每一行的數(shù)據(jù)都要執(zhí)行FOR循環(huán)的命令。其中 最關(guān)鍵的命令為DWT命令(離散小波)，這個命令把每一行波動 進(jìn)行小波分析并且提出其中的特色。其中我們使用的是其中HAAR小波，BOUNDARY設(shè)置為PERIODIC!就是循環(huán)的。RBIND命令只是把所有的特色困 綁在一起存儲于FEATUR度景中。最后一句的命令是將

7、生成的FEATURE陣轉(zhuǎn)化為DATA.FRAMEE寸象可做下面的使用。這里是上面一部分代碼運行的貼圖pClassId F exczaccing EBT coeffxciantJvxth Haar filter) librais-(vavelccs)pckag waveletsrw&g bJLlz under veiBicr; 2.14.2faatur*nydsra tnYdete for ii in 1 :nrowini7da3Ai l i tWL- rncarhaar, uounddiy=,feioai3)*recLUie - rblnlffaLar,f edtuze 1# set

8、class labels into categorical values classId - c(rep(T,100), rep(2,100),rep(3,100),rep(4,100), rep(5,100), rep(6,100) wtSc - data.frame(cbind(classId, feature)我們數(shù)據(jù)中有6種波，上面的命令是將波的類型(1-6)和波的特色捆綁在一起 生成新的DATA.FRAMEE寸象稱作wtSc。# build a decision tree with ctree() in package party library(party)ct (aux(cla3

9、sld=pla58ld) / nrow(wtSc)1371667plot(ct, ip_args=list(pval=FALSE), ep_args=list(digits=0)PLOT語句用來輸出決策樹模型的插圖回到這里從新做了心電圖。數(shù)據(jù)中1-100行都是1型波動。通過我們的數(shù)學(xué)模 型，還可以知道這個病人的心電圖是1型的嗎？通過數(shù)學(xué)模型來判定病人的心 電圖型531比方說，數(shù)據(jù)第一行代表的病人,或者是一個與這個病人狀況非常相似的病人,號被稱作預(yù)測。下面就是進(jìn)行預(yù)測的代碼，同樣需要將心電圖的特色提 取出來，然后放回到數(shù)學(xué)模型中來判定型號。#make predictionsf-NULL初始化F變

10、景a - t(mydata1,) 1號病人的心電圖存入a中，代表這個人又回來了。wt - dwt(a, filter=haar”, boundary=periodic)輸入提取特色f - rbind(f, unlist(c(wtW,wtVwtlevel)特色放入f中f tmaJce predictions f a wt HKLU旦(W,(G (Wt&W# Wt&1V1 J ) ) f predict(ct,f)：1 1Levels: 1 2 3 4 S 6預(yù)測結(jié)果是1 1是說明我們的數(shù)學(xué)模型已經(jīng)認(rèn)出這個病人的心電圖是1型的#/datab

11、ases/synthetic_control/synthetic_control.data.html#http:/ extracting DWT coefficients (with Haar filter)library(wavelets)feature-NULLmydata -read.table(/ml/databases/synthetic_control/synthetic_control.data”,header=F,sep=)#mydata - read.table(C:/Users/User/Desktop/RLanguage/

12、Wavelet/synthetic_control.data,header=F, sep=)for (i in 1:nrow(mydata) a - t(mydatai,)wt - dwt(a, filter=haar, boundary=periodic)feature - rbind(feature, unlist(c(wtW,wtVwtlevel)feature - as.data.frame(feature)# set class labels into categorical valuesclassId - c(rep(1,100), rep(2,100), rep(3,100),r

13、ep(4,100), rep(5,100), rep(6”,100)wtSc - data.frame(cbind(classId, feature)# build a decision tree with ctree() in package partylibrary(party)ct - ctree(classId ., data=wtSc,controls = ctree_control(minsplit=30, minbucket=10, maxdepth=5)pClassId - predict(ct)# check predicted classes against original class labelstable(classId, pClassId)# accuracy(sum(classId=pClassId) / nrow(wtSc