########　漁港フィルターの計算　（作成佐伯公康 H27年03月14日版）　##########

#★計算を始める前の注意事項★

#★１★計算を始める前に、ここから40行ほど下の、■１と書かれたすぐ下にある、bvalue、Da、N2T のパラメータの確認を行い、必要に応じて書き換えてください。

#★２★計算を始める前に、入出力ファイル名等の確認を行い、必要に応じて書き換えてください。具体的には次の通りです。
#　ここから50行ほど下の、■２Ａの直下の行で、計算に用いる加速度時刻歴のファイル名の指定を行います。このとき、フォルダ構成もあわせて指定します。
#　ここから50行ほど下の、■２Ｂの直下の行で、計算で得られた照査用震度などを出力するファイル名の指定を行います。このとき、フォルダ構成もあわせて指定します。
#　そのため、■２Ａの直下の行と、■２Ｂの直下の行について、お使いのＰＣのフォルダ構成に応じて書き換えてから使ってください。また、ファイル名は計算のたびごとに必要に応じて書き直して使ってください。

#★３★加速度時刻歴は次のようなデータとして、ファイルとして用意してください。
#・縦方向に（上から下に向けて）、一連の加速度時刻歴がずらりと並ぶデータを、テキスト形式のファイルとして用意してください。
#・データ数は2の累乗としてください。
#・上記の一列の加速度時刻歴以外のデータは、一切記さないでください。
#・たとえば、エクセルシートの、数字の入った縦一列のセルをコピーしてメモ帳に張り付け、txt形式で保存すると、上記のようなテキスト形式ファイルができます。
#・加速度の単位は、cm/s2としてください。
#・データの形式（小数形式、指数形式等）によってはRに適切に読み込めない場合があります。適切に読み込めない場合はデータの形式の変更をご検討ください。

#★以上のことを行いましたら、このプログラムの最初から最後まですべてを、Rに貼り付けてください。

#★なお、本プログラムを少し書き換えると、次のことも実行できます。くわしくは佐伯へお尋ねください。
#・加速度時刻歴を、ファイルからではなく、計算途中にエクセルからコピペで読み込む
#・フィルターを掛けた時系列波形を出力する
#　
#★注意事項　ここまで★



rm(list=ls(all=TRUE))   #←計算を始める前に、念のため、Rのなかに残っているすべての変数を消去するコマンド

#■１　次の行を書き換えて使う
bvalue <-  0.77 # b値（別途計算して、記入してください）
Da     <-  10   # 許容される重力式係船岸天端における変形量（cm）（設定した値を記入してください）
N2T <- 16384  # データ個数（加速度時刻歴のデータの長さを記入してください。つまり、後述する加速度時刻歴ファイルにおける、縦方向に並んだデータの数です。データ個数は、２の累乗にしてください）


ND <- 1    # 計算を行う波の数 （注）今のところ1以外は使えません(プログラム最後のパラメータ算定が、2以上に対応できない）
DT <- 0.01   # サンプリング周期（秒）（注）0.01以外にしますと、照査用震度の計算式を適用できなくなります。



#■２Ａ　次の行を書き換えて使う
input_f  <- ("C:/Users/Kimiy/Documents/R_calcdata/w01otsu_1.txt")
#■２Ｂ　次の行を書き換えて使う
output_f <- ("C:/Users/Kimiy/Documents/R_calcdata/fltdpara01_1.csv")


###　setwd("C:/Users/Kimiy/Documents/R output receptacle")


NT <- N2T/2+1  # 16384だと8193


DDY <- read.table(input_f, header=F)
## DDY <- read.table("clipboard",header=F)          #●●●#


#空の複素数ベクトルの作成
Sf       <- complex(NT)
fltdspec <- complex(NT)
cplxC    <- complex(N2T)
af       <- complex(NT)


#周波数刻みの作成
f1<- 1/N2T/DT  # 16384 と 0.01 なら 0.0061
kizami  <-seq(0,(1/DT/2), by=f1)  # 大崎p.39

#漁港フィルターの、フラットと曲線部の境界周波数の位置
Nfb <- round (1.2 / f1)



#時間刻みの作成
timesec <-seq(0,(N2T-1)*DT, by=DT)


Outfltd   <- data.frame(timesec)   #結果時系列を収めるデータフレームの初期型


#★繰り返し演算　ここから
for (ir in 1:ND){

#加速度時系列のベクトル化
DDYF<- numeric(N2T)                # DDYFの初期型（0を並べたベクトル）
DDYF[1:N2T]  <- DDY[, ir]          # DDYFへDDYのうち1列を代入

#フーリエ変換 ( FFT )
fftrslt <- fft  (DDYF)

#漁港フィルター作成と掛け算 Nfbを境界として低周波数帯と高周波数帯で別の式を使う。


for (ik in 1:NT)
   {
if ( Nfb >= ik ){
af[ik] <- bvalue
 } else {
gf <- 0.099*((ik-1)*f1 - 1.2)
af[ik] <- bvalue/(1 - gf^2 + 18.5*gf*(0+1i) )
}
fltdspec[ik] <- fftrslt[ik] * af[ik]
}



#  for (ik in 1:NT)
#   {
#  ff <- kizami[ik]
#  Sf[ik] <- 1 / (1-(ff/7.14)^2 + 2* (ff/7.14)*(0+1i))
#  fltdspec[ik] <- fftrslt[ik] * Sf[ik]
#   }

#漁港フィルターを1とした場合の演算（プログラミングの作業用）
# for (ik in 1:NT)
#   {
#   Sf[ik] <- 1
#   fltdspec[ik] <- fftrslt[ik] * Sf[ik]
#   }

#鏡像の作成
#cplxCの前半分の値を埋める
cplxC[1:NT]  <-  fltdspec  #1から8193

#cplxCの後半分の値を埋める
#繰り返しの冒頭
for(nhf in 2:(NT-1))       #2から8192
  {
 cplxC [(NT*2-nhf)] <- Conj(cplxC [nhf])    #16384から8194
  }

resulttsris  <- Re ( fft(cplxC, inverse = TRUE) / length(cplxC) )

Outfltd   <- data.frame(Outfltd,resulttsris)

}
#★繰り返し演算　ここまで


SSS <- sqrt(sum(resulttsris * resulttsris))
alphaf  <- max(abs(max(resulttsris)),abs(min(resulttsris)))
ppp <- 0.36*log(SSS/alphaf)-0.29
alphac  <- ppp*alphaf

kh  <- 1.78*(Da/10)^(-0.55)*alphac/980+0.04


#計算で求まったパラメータの画面表示

#↓フィルターを乗じた後の最大加速度 αf
alphaf
#↓継続時間による低減率p
ppp
#↓継続時間による補正をした後の最大加速度 αc
alphac
#↓照査用震度kh
kh

#ファイルへの書き出し
outparas <- data.frame( parameters=c(input_f,output_f,"b_value","Da","alpha_f","p","alpha_c","kh"), values=c( 9999,9999,bvalue,Da,alphaf,ppp,alphac,kh) )
write.csv(outparas, file=output_f, row.names=F)
# ここで、9999とは空欄を埋めるためのダミーの数値です。


## write.table(Outfltd ,file="outfltd.txt" ) #←これは、フィルターを掛けた時系列波形を出力させるコマンド。現状では使用不可。使用するには本プログラムの書き換えを要します。





#▼ここから下は、グラフ作成コマンド。必要に応じてX,Y軸範囲などを書き換えてご使用ください。

par ( bg = "grey 95",  mar = c ( 4,4,1,1 ),  mfrow = c ( 2, 2 ))


plot(timesec[],DDY[,1],xlim=c(0,170),ylim=c(-300,300), type="l",
　　　xlab="Time(sec)",ylab="Accerelation(cm/s2)",col="black")     

mtext("上（黒）：フィルターを掛ける前の加速度時系列。下（赤）：フィルターを掛けた後の加速度時系列。右（青）：フィルター関数。", side = 3, line = 0, at = NA)

realaf <- Re(af)
plot(kizami[],realaf[],xlim=c(0.1,10),ylim=c(0,1.6), type="l",
　　　xlab="Frequency(Hz)",ylab="a(f)",col="blue",log="x" )     


plot(timesec[],resulttsris[],xlim=c(0,170),ylim=c(-300,300), type="l",
　　　xlab="Time(sec)",ylab="Accerelation(cm/s2)",col="red")