Sub 震度階() Dim xr() As Double, xi() As Double Dim yr() As Double, yi() As Double Dim zr() As Double, zi() As Double Dim temp As Double, a As Double, s As Double, dt As Double Dim n As Long, nn As Long, nfold As Long, i As Long Dim ind As Integer n = [a2] dt = [a4] temp = Log(n) / Log(2) nn = 2 ^ CInt(temp) nfold = nn / 2 + 1 ReDim xr(nn), xi(nn), yr(nn), yi(nn), zr(nn), zi(nn) 'データの読み込み For i = 1 To n xr(i) = Cells(i + 1, 2): xi(i) = 0# yr(i) = Cells(i + 1, 3): yi(i) = 0# zr(i) = Cells(i + 1, 4): zi(i) = 0# Next i For i = n + 1 To nn xr(i) = 0#: xi(i) = 0# yr(i) = 0#: yi(i) = 0# zr(i) = 0#: zi(i) = 0# Next i 'FFTで周波数領域に ind = -1 Call fast(xr, xi, nn, ind) Call fast(yr, yi, nn, ind) Call fast(zr, zi, nn, ind) 'フィルターをかける Call filter(xr, xi, dt, nn) Call filter(yr, yi, dt, nn) Call filter(zr, zi, dt, nn) 'FFTで時間領域に戻す ind = 1 Call fast(xr, xi, nn, ind) Call fast(yr, yi, nn, ind) Call fast(zr, zi, nn, ind) For i = 1 To n xr(i) = xr(i) / CDbl(nn) yr(i) = yr(i) / CDbl(nn) zr(i) = zr(i) / CDbl(nn) Next i '3成分をベクトル合成してE列に書き出す For i = 1 To n Cells(i + 1, 5) = Sqr(xr(i) ^ 2 + yr(i) ^ 2 + zr(i) ^ 2) Next i '大きい順に並び替え Range(Cells(2, 5), Cells(n + 1, 5)).Sort key1:=Range("e1"), order1:=xlDescending '0.3/dt番目の値aを求める i = Int(0.3 / dt) a = Cells(i + 1, 5) 's = 2 log a + 0.94 から計測震度sを求める(小数第3位を四捨五入し,第2位を切り捨て) s = 2 * Log(a) / Log(10) + 0.94 s = Int(Int(100# * s + 0.5) / 10#) / 10# Cells(2, 6) = s '震度階級 If s < 0.5 Then Cells(3, 6) = "震度0" ElseIf s < 1.5 Then Cells(3, 6) = "震度1" ElseIf s < 2.5 Then Cells(3, 6) = "震度2" ElseIf s < 3.5 Then Cells(3, 6) = "震度3" ElseIf s < 4.5 Then Cells(3, 6) = "震度4" ElseIf s < 5 Then Cells(3, 6) = "震度5弱" ElseIf s < 5.5 Then Cells(3, 6) = "震度5強" ElseIf s < 6 Then Cells(3, 6) = "震度6弱" ElseIf s < 6.5 Then Cells(3, 6) = "震度6強" Else Cells(3, 6) = "震度7" End If End Sub '------------フーリエ変換 Sub fast(xr() As Double, xi() As Double, nn As Long, ind As Integer) ' 参考文献:大崎順彦「新・地震動のスペクトル解析入門」鹿島出版会 Dim tr As Double, ti As Double Dim theta As Double Dim i As Long, j As Long, k As Long Dim m As Long, kmax As Long, istep As Long Dim PI As Double PI = Application.WorksheetFunction.PI() j = 1 For i = 1 To nn If i < j Then tr = xr(j): ti = xi(j) xr(j) = xr(i): xi(j) = xi(i) xr(i) = tr: xi(i) = ti End If m = nn / 2 Do While (j > m) j = j - m m = m / 2 If (m < 2) Then Exit Do Loop j = j + m Next i kmax = 1 Do While kmax < nn istep = kmax * 2 For k = 1 To kmax theta = PI * CDbl(ind * (k - 1)) / CDbl(kmax) For i = k To nn Step istep j = i + kmax tr = xr(j) * Cos(theta) - xi(j) * Sin(theta) ti = xr(j) * Sin(theta) + xi(j) * Cos(theta) xr(j) = xr(i) - tr: xi(j) = xi(i) - ti xr(i) = xr(i) + tr: xi(i) = xi(i) + ti Next i Next k kmax = istep Loop End Sub '---------フィルター Sub filter(xr() As Double, xi() As Double, dt As Double, nn As Long) Dim nfold As Long, i As Long Dim f As Double, y As Double Dim f1 As Double, f2 As Double, f3 As Double nfold = nn / 2 + 1 xr(1) = 0#: xi(1) = 0# For i = 2 To nfold - 1 f = CDbl(i - 1) / CDbl(nn) / dt y = f / 10# f1 = Sqr(1# / f) f2 = 1# / Sqr(1# + 0.694 * y ^ 2 + 0.241 * y ^ 4 + 0.0557 * y ^ 6 _ + 0.009664 * y ^ 8 + 0.00134 * y ^ 10 + 0.000155 * y ^ 12) f3 = Sqr(1# - Exp(-(2# * f) ^ 3)) xr(i) = f1 * f2 * f3 * xr(i) xi(i) = f1 * f2 * f3 * xi(i) xr(nn - i + 2) = xr(i) xi(nn - i + 2) = -xi(i) Next i f = CDbl(nfold - 1) / CDbl(nn) / dt y = f / 10# f1 = Sqr(1# / f) f2 = 1# / Sqr(1# + 0.694 * y ^ 2 + 0.241 * y ^ 4 + 0.0557 * y ^ 6 _ + 0.009664 * y ^ 8 + 0.00134 * y ^ 10 + 0.000155 * y ^ 12) f3 = Sqr(1# - Exp(-(2# * f) ^ 3)) xr(nfold) = f1 * f2 * f3 * xr(nfold) xi(nfold) = f1 * f2 * f3 * xi(nfold) End Sub