利用numpy計算傅裡葉譜和地震動傅裡葉振幅譜

引言

傅裡葉譜和地震動中的傅裡葉振幅譜是有些許差别的。主要表現在幅值的調整上。

參考資料：《地震動的譜分析入門》大崎順彥

傅裡葉譜

一般的傅裡葉譜的幅值利用fft變換後乘以2再除以N進行調整。

import numpy as np
import math
import matplotlib.pyplot as plt

sampRat = 100
T = 6

t = np.linspace(0, T, T*sampRat, endpoint=False)
y = 5*np.sin(2*math.pi*t*5)+10*np.sin(2*math.pi*t*25)

plt.figure()
plt.title('時程曲線')
plt.plot(t, y)
plt.xlabel('時間')
plt.ylabel('幅度')
plt.show()

f = np.linspace(0, sampRat, T*sampRat, endpoint=False)
ff = np.fft.fft(y)
ff = np.abs(ff)
ff = ff*2/sampRat/T

plt.figure()
plt.title('傅裡葉譜')
plt.plot(f, ff)
plt.show()
           

地震動傅裡葉振幅譜

地震動傅裡葉振幅譜利用fft變換後除以采樣頻率即可。

# -*- coding: utf-8 -*-

import numpy as np
import math
import matplotlib.pyplot as plt

sampRat = 100
T = 6

t = np.linspace(0, T, T*sampRat, endpoint=False)
y = 5*np.sin(2*math.pi*t*5)+10*np.sin(2*math.pi*t*25)

plt.figure()
plt.title('時程曲線')
plt.plot(t, y)
plt.xlabel('時間')
plt.ylabel('幅度')
plt.show()

f = np.linspace(0, sampRat, T*sampRat, endpoint=False)
ff = np.fft.fft(y)
ff = np.abs(ff)
ff = ff/sampRat

plt.figure()
plt.title('傅裡葉振幅譜')
plt.plot(f, ff)
plt.show()