R

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

Rcpp パッケージ

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

Rcpp 入門

OSはzorin linuxです。

Hanning Window の関数を定義(Rcppを利用)

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
library(Rcpp)
cppFunction('
            NumericVector hanningRcpp(List x) 
{          
             NumericVector acc=x["acc"];
            int nhan=x["nhan"];
            int nd = acc.size();
            NumericVector ret(nd);
//        NumericVector temp(nd);
//        for(int i=0; i< nd; i++){
//               temp[i]=acc[i]; 
//        }
            NumericVector temp=clone(acc);
            for(int ihan=0; ihan < nhan; ihan++){
            ret[0]=0.5*(temp[0]+temp[1]); 
            ret[nd-1]=0.25*temp[nd-2]+0.5*temp[nd-1]+0.25*temp[0]; 
            for(int i=1; i < nd-1; i++){
            ret[i]=0.25*temp[i-1]+0.5*temp[i]+0.25*temp[i+1]; 
            }
//            for(int i=0; i < nd; i++){
//            temp[i]=ret[i]; 
//           }
            temp=clone(ret);
            }
            return ret;
}
            ')

データを取り込み、メイン処理
地震動のスペクトル1と異なるのは１箇所（処理に使うHanning Window の関数が異なる）だけ

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
#入力
#加速度　境港市東本町
seibu1<-read.csv("http://www.data.jma.go.jp/svd/eqev/data/kyoshin/jishin/001006_tottori-seibu/dat/AA06E9E1.csv",skip=6,header=T)
#サンプリング周期
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)) 
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)}
}

グラフは「地震動のスペクトル1」と同じなので省略

計算の実行時間は約10分の１に減りました。（環境による）