R

地震動のスペクトル1(Hanningウィンドウ)

quantmod , RSEIS パッケージ

（データ）気象庁
強震波形（平成12年(2000年)鳥取県西部地震）の
境港市東本町
（参考）
フーリエスペクトル解析
Octaveによるフーリエ解析
WANtaroHP (R tips)の「フーリエスペクトル」（SWIN関数:Parzenウィンドウ）

OSはzorin linuxです。

Hanning Window の関数を定義

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
hanningR<-function(pd,nhan){
y<-NULL;yp1<-NULL;ym1<-NULL;temp<-NULL
y<-pd
yp1<-pd[2:length(y)]
yp1[length(y)]<-pd[1]
ym1[1]<-pd[2]
ym1[2:length(y)]<-pd[1:(length(y)-1)]
for (i in 1:nhan) {
temp<-0.25*ym1+0.5*y+0.25*yp1
y<-temp
yp1<-temp[2:length(temp)]
yp1[length(temp)]<-temp[1]
ym1[1]<-temp[2]
ym1[2:length(temp)]<-temp[1:(length(temp)-1)]
}
return (temp)
}

データを取り込み波形をplot

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#加速度　境港市東本町
seibu1<-read.csv("http://www.data.jma.go.jp/svd/eqev/data/kyoshin/jishin/001006_tottori-seibu/dat/AA06E9E1.csv",skip=6,header=T)
#波形を見る(ts class)
seibu1.ts<-ts(seibu1,start=0,freq=100)
head(seibu1.ts)
max(seibu1);min(seibu1)
#png("hakei01.png",width=1000,height=800)
par(mfrow=c(3,1))
plot(seibu1.ts[,1],ylim=c(-750,750),ylab="NS")
title("波形（加速度）")
plot(seibu1.ts[,2],ylim=c(-750,750),ylab="EW")
plot(seibu1.ts[,3],ylim=c(-750,750),ylab="UD")
par(mfrow=c(1,1))
#dev.off()

メイン

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#サンプリング周期
dt     <- 0.01
#windowのバンド幅
bd<-1
for (direction in 1:3){
if (direction==1) {d<-seibu1$NS}
if (direction==2) {d<-seibu1$EW}
if (direction==3) {d<-seibu1$UD}
    k  <- length(d)
#入力データ数が２の乗数になるように調整
    nn <- 2
    while( nn < k) nn <- nn*2
    x      <- numeric(nn)
    x[1:k] <- d[1:k]
    z      <- fft(x)
    z1     <- Re(z)/nn
    z2     <- Im(z)/nn
    nfold  <- nn/2+1
    df     <- 1.0/dt/nn
    fs1    <- numeric(nfold)
    fs2    <- numeric(nfold)
    frq    <- numeric(nfold)
    fs1    <- dt*nn*sqrt(z1[1:nfold]^2+z2[1:nfold]^2)
    ii     <- seq(0,nfold-1,by=1)
    frq    <- ii*df
#ウィンドウのバンド幅 bdのときの繰り返し回数nを計算 
hann<-(bd*(3*(length(x)*dt)/8))^2
n=ceiling(hann)
band<-round(8/3*sqrt(n)*(1/(length(x)*dt)),3)
fs2<-hanningR(fs1,n)
#fs2<-hanningRcpp(list(acc=fs1,nhan=n)) 
#options(scipen=10)
if (direction==1) {resNS<-data.frame(frq,fs1,fs2)}
if (direction==2) {resEW<-data.frame(frq,fs1,fs2)}
if (direction==3) {resUD<-data.frame(frq,fs1,fs2)}
}

plot1

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
#png("spec01.png",width=1000,height=800)
#gridx=c(seq(2*dt,0.09,0.01),seq(0.1,0.9,0.1),seq(1,10,1))
gridx=c(seq(0,20,1))
gridy=c(seq(0.001,0.009,0.001),seq(0.01,0.09,0.01),seq(0.1,0.9,0.1),seq(1,9,1),seq(10,90,10),seq(100,900,100),seq(1000,10000,1000))
par(mfrow=c(1,3))
plot(resNS$frq,resNS$fs1,type="l",log="y",las=1,col="red",xlab="Frequency (Hz)",ylab="",yaxt="n",
panel.first = abline(v=gridx,h=gridy,col="lightgray",lty ="dotted"),xlim=c(0,20),ylim=c(1e-2,7e+2))
axis(2,c(0.00001,0.0001,0.001,0.01,0.1,1,10,100,1000,10000),
c("0.00001","0.0001","0.001","0.01","0.1","1","10","100","1000","10000"),las=1) 
lines(resNS$frq,resNS$fs2,col="black",lwd=2)
title("Fourier spectrum (gal * sec) NS")
legend(
    "bottomleft",
    c("Original",paste("Hanning Window :band=",band,"Hz : n=",n)),
    col = c("red","black"),
    lty = 1,
    lwd = 2,
    bg  = "white"
)
plot(resEW$frq,resEW$fs1,type="l",log="y",las=1,col="red",xlab="Frequency (Hz)",ylab="",yaxt="n",
panel.first = abline(v=gridx,h=gridy,col="lightgray",lty ="dotted"),xlim=c(0,20),ylim=c(1e-2,7e+2))
axis(2,c(0.00001,0.0001,0.001,0.01,0.1,1,10,100,1000,10000),
c("0.00001","0.0001","0.001","0.01","0.1","1","10","100","1000","10000"),las=1) 
lines(resEW$frq,resEW$fs2,col="black",lwd=2)
title("Fourier spectrum (gal * sec) EW")
legend(
    "bottomleft",
    c("Original",paste("Hanning Window :band=",band,"Hz : n=",n)),
    col = c("red","black"),
    lty = 1,
    lwd = 2,
    bg  = "white"
)
plot(resUD$frq,resUD$fs1,type="l",log="y",las=1,col="red",xlab="Frequency (Hz)",ylab="",yaxt="n",
panel.first = abline(v=gridx,h=gridy,col="lightgray",lty ="dotted"),xlim=c(0,20),ylim=c(1e-2,7e+2))
axis(2,c(0.00001,0.0001,0.001,0.01,0.1,1,10,100,1000,10000),
c("0.00001","0.0001","0.001","0.01","0.1","1","10","100","1000","10000"),las=1) 
lines(resUD$frq,resUD$fs2,col="black",lwd=2)
title("Fourier spectrum (gal * sec) UD")
legend(
    "bottomleft",
    c("Original",paste("Hanning Window :band=",band,"Hz : n=",n)),
    col = c("red","black"),
    lty = 1,
    lwd = 2,
    bg  = "white"
)
par(mfrow=c(1,1))
#dev.off()

plot2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#png("spec02.png",width=1000,height=800)
gridx=c(seq(2*dt,0.09,0.01),seq(0.1,0.9,0.1),seq(1,10,1))
#gridx=c(seq(0,20,1))
gridy=c(seq(0.001,0.009,0.001),seq(0.01,0.09,0.01),seq(0.1,0.9,0.1),seq(1,9,1),seq(10,90,10),seq(100,900,100),seq(1000,10000,1000))
plot(1/resEW$frq[-1],resEW$fs2[-1],type="l",log="xy",las=1,col="red",xlab="Period(sec)",ylab="",yaxt="n",
panel.first = abline(v=gridx,h=gridy,col="lightgray",lty ="dotted"),xlim=c(0.027,10),ylim=c(0.1,500))
axis(2,c(0.00001,0.0001,0.001,0.01,0.1,1,10,100,1000,10000),
c("0.00001","0.0001","0.001","0.01","0.1","1","10","100","1000","10000"),las=1) 
lines(1/resNS$frq[-1],resNS$fs2[-1],col="blue")
lines(1/resUD$frq[-1],resUD$fs2[-1],col="green")
title(paste("Fourier spectrum (gal * sec) [band=",band,"Hz]"))
legend(
    "bottomright",
    c("NS","EW","EW"),
    col = c("blue","red","green"),
    lty = c(1,1,1),
    lwd = c(2,2,2),
    bg  = "white"
)
#dev.off()

plot3

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
#png("spec03.png",width=1000,height=800)
gridx=c(seq(0,20,1))
gridy=c(seq(0.001,0.009,0.001),seq(0.01,0.09,0.01),seq(0.1,0.9,0.1),seq(1,9,1),seq(10,90,10),seq(100,900,100),seq(1000,10000,1000))
plot(resNS$frq,resNS$fs1,type="l",log="y",las=1,col="red",xlab="Frequency (Hz)",ylab="",yaxt="n",
panel.first = abline(v=gridx,h=gridy,col="lightgray",lty ="dotted"),xlim=c(0,20),ylim=c(1e-2,7e+2))
axis(2,c(0.00001,0.0001,0.001,0.01,0.1,1,10,100,1000,10000),
c("0.00001","0.0001","0.001","0.01","0.1","1","10","100","1000","10000"),las=1) 
lines(resNS$frq,resNS$fs2,col="black",lwd=2)
title("Fourier spectrum (gal * sec) NS")
legend(
    "bottomleft",
    c("Original",paste("Hanning Window :band=",band,"Hz : n=",n)),
    col = c("red","black"),
    lty = 1,
    lwd = 2,
    bg  = "white"
)
#ピーク検出
#quantmodパッケージのfindPeaks関数
library(quantmod)
peak<-data.frame(frq=resNS$frq,fs2=resNS$fs2)
p=findPeaks(subset(peak, subset=frq<=21)$fs2)
points(peak$frq[p],peak$fs2[p], col=rgb(0,0,1,0.7),pch=20,cex=3)
#dev.off()
library(knitr)
kable(data.frame(Frequency=peak$frq[p],Peak=peak$fs2[p]))

Frequency Peak
0.6500244 113.352854
2.3406982 108.850808
3.4667969 69.430107
4.9743652 55.929940
8.6975098 10.572927
10.9344482 8.278958
12.4053955 3.793472
13.5284424 4.100282
15.3930664 2.175454
16.9860840 2.018440
18.2647705 1.555040
19.7082520 1.035607
20.8465576 1.068462

速度波形・変位波形の求め方によると
「
速度波形はカットオフ5秒の３次のバターワースフィルターを用いて積分を行っている。」とのこと

これであっているの？？？

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
library(RSEIS)
xNS <- butfilt(seibu1.ts[,1],fl=0.2,fh=50,deltat=0.01,"HP","BU",npoles=3)
velNS<-0.01*cumsum(xNS) ; velNS.ts<-ts(velNS,start=0,freq=100)
xEW <- butfilt(seibu1.ts[,2],fl=0.2,fh=50,deltat=0.01,"HP","BU",npoles=3)
velEW<-0.01*cumsum(xEW) ; velEW.ts<-ts(velEW,start=0,freq=100)
xUD <- butfilt(seibu1.ts[,3],fl=0.2,fh=50,deltat=0.01,"HP","BU",npoles=3)
velUD<-0.01*cumsum(xUD) ; velUD.ts<-ts(velUD,start=0,freq=100)
#png("hakei02.png",width=1000,height=800)
par(mfrow=c(3,1))
plot(velNS.ts,ylab="NS",ylim=c(-40,40))
title("波形（速度）")
plot(velEW.ts,ylab="EW")
plot(velUD.ts,ylab="UD",ylim=c(-40,40))
par(mfrow=c(3,1))
#dev.off()