成果展示

第8章時(shí)間序列的譜分析WUSTY.M.Ding

本章結(jié)構(gòu)頻域1.時(shí)間序列的周期性2.譜密度3.譜估計(jì)4.線性濾波器5.問(wèn)題的提出人類(lèi)視覺(jué)所感受到的是在空間域和時(shí)間域的信號(hào)。但是，往往許多問(wèn)題在頻域中討論時(shí)，有其非常方便分析的一面。例如，空間位置上的變化不改變信號(hào)的頻域特性。

如何理解空間域、時(shí)間域與頻域？問(wèn)題的提出：如何理解頻域？

音樂(lè)最普遍的理解：隨著時(shí)間變化的震動(dòng)

音樂(lè)更直觀的理解空域頻域/p/19763358

同一物體的不同視角理解問(wèn)題的提出

時(shí)域：以時(shí)間作為自變量(例子：波形圖)

頻域：以頻率作為自變量(例子：頻譜圖)

空間域：以空間坐標(biāo)作為自變量(例子：數(shù)字圖像)

行（i）列（j）矩陣A(i,j)圖像矩陣坐標(biāo)系右圖為同一信號(hào)，相比之下，頻域分析該信號(hào)更加方便問(wèn)題的提出

時(shí)域分析只能反映信號(hào)的幅值隨時(shí)間的變化情況

頻域分析反映了信號(hào)不同頻率組成及其分量成分大小，能夠提供比時(shí)域信號(hào)波形更直觀，豐富的信息。

圖例：受噪聲干擾的多頻率成分信號(hào)

基于頻域的基本組成單元-正弦波，一個(gè)矩形波在頻域中的表示如下圖所示：

w傅立葉變換理解：歐拉公式

歐拉公式關(guān)鍵作用：將正弦波統(tǒng)一成了簡(jiǎn)單的指數(shù)形式

問(wèn)題：正弦波形式復(fù)雜，周期函數(shù)，不易數(shù)學(xué)與計(jì)算機(jī)處理

解決方法：歐拉公式傅立葉變換理解：一維FT及其反變換

連續(xù)函數(shù)f(x)的傅立葉變換F(u):

傅立葉變換F(u)的反變換：

將歐拉公式帶入可得離散傅立葉公式：（8.1）（8.2）（8.3）傅立葉變換理解：一維FT及其反變換

F(u)由f(x)與對(duì)應(yīng)頻率的正弦和余弦乘積和組成

u值決定了變換的頻率成份，因此，F(xiàn)(u)覆蓋的域稱為頻率域，其中每一項(xiàng)都被稱為FT的頻率分量

頻率域及頻率分量與f(x)的時(shí)間域和時(shí)間成分相對(duì)應(yīng)（8.4）

本章結(jié)構(gòu)頻域1.時(shí)間序列的周期性2.譜密度3.譜估計(jì)4.線性濾波器5.

時(shí)間序列的周期性

時(shí)間序列的周期性

時(shí)間序列的周期性

時(shí)間序列的周期性一個(gè)例子是，在電影中，攝像機(jī)對(duì)移動(dòng)汽車(chē)上的旋轉(zhuǎn)車(chē)輪進(jìn)行采樣，其中車(chē)輪似乎以不同的速度旋轉(zhuǎn)。電影以每秒24幀的速度錄制。如果相機(jī)拍攝的是一個(gè)以24周/秒(或24赫茲)的速度旋轉(zhuǎn)的輪子，這個(gè)輪子看起來(lái)是靜止的。這個(gè)鏈接可以了解更多相關(guān)信息/watch?v=6XwgbHjRo30

時(shí)間序列的周期性

時(shí)間序列的周期性

時(shí)間序列的周期性

時(shí)間序列的周期性

時(shí)間序列的周期性

時(shí)間序列的周期性

時(shí)間序列的周期性

P(40/100)=P(60/100)=85,其他點(diǎn)處P(j/n)=0。這些值正好是例子中產(chǎn)生的分量振幅的平方值。我們從前面的案例中模擬數(shù)據(jù)x,R代碼如下n=length(x)P=(4/n)*Mod(fft(x)/sqrt(n))^2Fr=0:99/100plot(Fr,P,type="o",xlab="frequency",ylab="scaledperiodogram")

時(shí)間序列的周期性

時(shí)間序列的周期性如果我們把例8.1中的數(shù)據(jù)xt看作是由不同強(qiáng)度(振幅)的三原色xt1、xt2、xt3組成的顏色(波形)，那么我們可以把周期圖看作是將顏色xt分解為其三原色(光譜)的棱鏡。因此有了頻譜分析這個(gè)術(shù)語(yǔ)。

例下圖的數(shù)據(jù)是在午夜連續(xù)600天拍攝的一顆恒星的級(jí)別。數(shù)據(jù)來(lái)源于E.T.惠特克和G.羅賓遜(1923,Blackie&Son,Ltd.)的經(jīng)典文本《CalculusofObservations,aTreatiseonNumericalMathematics》，圖中還顯示了小于.08頻率的周期圖;高于0.08頻率的周期圖基本上為零。請(qǐng)注意，29天周期和25天周期是數(shù)據(jù)中最顯著的周期成分。相應(yīng)的R代碼：

時(shí)間序列的周期性n=length(star)par(mfrow=c(2,1),mar=c(3,3,1,1),mgp=c(1.6,.6,0))tsplot(star,ylab="starmagnitude",xlab="day")Per=Mod(fft(star-mean(star)))^2/nFreq=(1:n-1)/nplot(Freq[1:50],Per[1:50],type='h',lwd=3,ylab="Periodogram",xlab="Frequency")u=which.max(Per[1:50])#22freq=21/600=.035cycles/dayuu=which.max(Per[1:50][-u])#25freq=25/600=.041cycles/day1/Freq[22];1/Freq[26]#period=days/cycletext(.05,7000,"24daycycle");text(.027,9000,"29daycycle")###anotherwaytofindthetwopeaksistoorderonPery=cbind(1:50,Freq[1:50],Per[1:50]);y[order(y[,3]),]

時(shí)間序列的周期性

本章結(jié)構(gòu)頻域1.時(shí)間序列的周期性2.譜密度3.譜估計(jì)4.線性濾波器5.8.3.1譜密度

8.3.1譜密度

8.3.1譜密度

8.3.2ARMA的譜密度MA(1)模型的譜密度

MA(1)模型的譜密度

圖8.8理論白噪聲頻譜(上)，一階移動(dòng)平均(中)，和二階自回歸過(guò)程(下)

R代碼實(shí)現(xiàn)：w=rnorm(550,0,1)#50extratoavoidstartupproblemsx=filter(w,filter=c(1,-.9),method="recursive")[-(1:50)]tsplot(x,main="autoregression")

本章結(jié)構(gòu)頻域1.時(shí)間序列的周期性2.譜密度3.譜估計(jì)4.線性濾波器5.

圖8.10圖8.9所示的AR(2)譜的一個(gè)小截面(接近峰值)

例8.8SOI和Recuitment的平均周期圖圖8.11SOI和Recruitment的周期圖，n=453(n0=480)，其中頻率軸以年為單位標(biāo)記。注意在ω=1處的共同峰值，或每年一個(gè)周期，以及ω=1/4附近的一些較大的值，或每四年一個(gè)周期圖8.12SOI和Recruitment系列的平均周期圖n=453,n0=480,L=9,df=17，在4年周期，ω=1周期/4年，年周期，ω=1周期/年，以及k=2,3時(shí)的一些諧波ω=k。圖8.12的R代碼：par(mfrow=c(2,1))(k=kernel("daniell",4))soi.ave=mvspec(soi,k,log="no");abline(v=c(.25,1,2,3),lty=2)rec.ave=mvspec(rec,k,log="no");abline(v=c(.25,1,2,3),lty=2)soi.ave$bandwidth#=0.225soi.ave$df

圖8.13

承接圖8.12，以log10尺度繪制的平均周期圖坐標(biāo)。右上角的顯示表示一般的95%置信區(qū)間例8.8諧波

諧波的R代碼為t=seq(0,1,by=1/200)amps=c(1,.5,.4,.3,.2,.1)x=matrix(0,201,6)for(jin1:6)x[,j]=amps[j]*sin(2*pi*t*2*j)x=ts(cbind(x,rowSums(x)),start=0,deltat=1/200)ts.plot(x,lty=c(1:6,1),lwd=c(rep(1,6),2),ylab="Sinusoids")names=c("Fundamental","2ndHarmonic","3rdHarmonic","4thHarmonic","5thHarmonic","6thHarmonic","FormedSignal")legend("topright",names,lty=c(1:6,1),lwd=c(rep(1,6),2))

圖8.14調(diào)整的Daniell內(nèi)核權(quán)值對(duì)應(yīng)的R代碼如下：k=kernel("modified.daniell",c(3,3))soi.smo=mvspec(soi,k,taper=.1,log="no")#ataperisusedabline(v=c(1/4,1),lty="dotted")##Repeatabovelineswithrecreplacingsoiinline3df=soi.smo$df#df=17.42618soi.smo$bandwidth#Bw=0.2308103=12*9.232/480圖8.15SOI和recruit系列的平滑(錐形)譜估計(jì)8.4.2參數(shù)譜估計(jì)

例SOI的自回歸譜估計(jì)量考慮獲得與SOI系列的圖所示的非參數(shù)估計(jì)量相比較的結(jié)果。對(duì)p=1、2、3時(shí)的高階AR(p)模型進(jìn)行逐次擬合，在p=15時(shí)得到最小的BIC和最小的AIC，如圖8.16所示。圖8.16模型選擇準(zhǔn)則AIC和BIC是擬合SOI系列自回歸模型的p階函數(shù)spaic=spec.ar(soi,log="no")#minAICspecabline(v=frequency(soi)*1/52,lty="dotted")#ElNinoCycle(soi.ar=ar(soi,order.max=30))#estimatesandAICsdev.new()plot(1:30,soi.ar$aic[-1],type="o")#plotAICsn=length(soi)c()->AIC->AICc->BICfor(kin1:30){sigma2=ar(soi,order=k,aic=FALSE)$var.predBIC[k]=log(sigma2)+(k*log(n)/n)AICc[k]=log(sigma2)+((n+k)/(n-k-2))AIC[k]=log(sigma2)+((n+2*k)/n)}IC=cbind(AIC,BIC+1)ts.plot(IC,type="o",xlab="p",ylab="AIC/BIC")

本章結(jié)構(gòu)頻域1.時(shí)間序列的周期性2.譜密度3.譜估計(jì)4.線性濾波器5.

圖8.17一階差分（上）核十二個(gè)月移動(dòng)平均（下）濾波后的平方頻率響應(yīng)函數(shù)

par(mfrow=c(3,1),mar=c(3,3,1,1),mgp=c(1.6,.6,0))tsplot(soi)#plotdatatsplot(diff(soi))#plotfirstdifferencek=kernel("modified.daniell",6)#filterweightstsplot(soif<-kernapply(soi,k))#plot12monthfilterdev.new