python绘制地震动傅里叶谱

最新推荐文章于 2023-12-23 14:29:33 发布

酒桶在你野区

最新推荐文章于 2023-12-23 14:29:33 发布

阅读量1k

点赞数 1

文章标签： python

本文链接：https://blog.csdn.net/s7777777777777/article/details/125956239

版权

先说结论

快速傅里叶变换结果为 $y_k$ ，我们需要的傅里叶谱振幅为 $F_k$ ，二者的关系如下
$F_k =\Delta t\cdot|y_k|$
需要注意： $F_0 = \Delta t/2\cdot|y_0|$ ， $F_{N/2} = \Delta t/2\cdot|y_{N/2}|$

python代码

import numpy as np
from scipy.fft import fft


def fourier(ag, dt):
    n = len(ag) - 1
    fw = fft(ag)
    i_half = n // 2
    fw_half = fw[:i_half + 1]  # 取单边频谱
    ampl = np.abs(fw_half) * dt  # 幅值
    ampl[0] = 0.  # 直流分量置零
    ampl[-1] = ampl[-1] / 2
    freq = np.linspace(0, 1 / (2 * dt), i_half + 1)
    return ampl, freq

原因

1 地震动的傅里叶变换与逆变换

地震动 $x_m$ 是离散值，傅里叶逆变换为
$x_m=\frac{A_{0}}{2}+\sum_{k=1}^{N / 2-1}\left[A_{k} \cos \frac{2 \pi k m}{N }+B_{k} \sin \frac{2 \pi k m}{N }\right]+\frac{A_{N / 2}}{2} \cos \frac{2 \pi(N / 2) m}{N }\tag{1}$
傅里叶变换为
$A_{k}=\frac{2}{N} \sum_{m=0}^{N-1} x_{m} \cos \frac{2 \pi k m}{N} \quad k=0,1,2, \cdots, \frac{N}{2}-1, \frac{N}{2}\tag{2}$
$B_{k}=\frac{2}{N} \sum_{m=0}^{N-1} x_{m} \sin \frac{2 \pi k m}{N} \quad k=1,2, \cdots, \frac{N}{2}-1 \tag{3}$
$A_k$ 为有限傅里叶余弦系数， $B_k$ 为有限傅里叶正弦系数， $A_k$ 比 $B_k$ 多两个。上述公式的推导过程不在此赘述。

2 复数表示的傅里叶变换和逆变换

根据欧拉公式 $\mathrm{e}^{\pm \mathrm{i} \theta}=\cos \theta \pm \mathrm{i} \sin \theta$ ，可以转化傅里叶逆变换式（1）为复数形式
$x_{m}=\frac{1}{2} \sum_{k=0}^{N-1}\left(A_{k}-\mathrm{i} B_{k}\right) \mathrm{e}^{\mathrm{i}(2 \pi k m / N)} \tag{4}$
可以得到复傅里叶系数和有限傅里叶系数的关系
$C_k=\frac{A_k-\mathrm{i}B_k}{2}\quad k=0,1,2, \cdots, N-1 \tag{5}$
相应的复数形式的傅里叶变换为
$C_{k}=\frac{1}{N} \sum_{m=0}^{N-1} x_{m} \mathrm{e}^{-\mathrm{i}(2 \pi k m / N)} \quad k=0,1,2, \cdots, N-1 \tag{6}$

3 python中的快速傅里叶变换

python中scipy库快速傅里叶变换的结果是下式中的 $y_k$
$y_k=\sum_{m=0}^{N-1} x_m e^{-i(2 \pi km/N)} \tag{7}$
对比式（7）和式（8），可以发现二者相差 $1/ N$

4 傅里叶谱振幅与相位

根据三角函数的关系式 $\cos \alpha+b \sin \alpha=\sqrt{a^{2}+b^{2}} \cos (\alpha+\phi)$ ，由式（1）得，地震动傅里叶谱的振幅 $X_k$ 为
$X_{k}=\sqrt{A_{k}^{2}+B_{k}^{2}}\quad k=1,2,\cdots,N/2-1 \tag{8}$
需要注意： $X_0=A_0/2$ 、 $X_{N/2}=A_{N/2}/2$

相位为
$\phi_{k}=\arctan \left(-\frac{B_{k}}{A_{k}}\right) \quad-\pi<\phi_{k}<\pi\tag{9}$
需要注意： $\phi_0=0$ 、 $\phi_{N/2}=0$

5 傅里叶积分

上面的所有处理初看起来是处理有限时间 $T$ 内的现象，实际上是处理周期为 $T$ 的无限次反复现象。
环状效应
而我们需要处理的地震动是在某一时刻开始，持续一段时间后就结束的不具有重复性的现象。因此，把处理对象扩展为定义于 $-\infty \leqslant t \leqslant \infty$ 的非周期函数 $x (t)$

$x (t)$ 的傅里叶变换（也即傅里叶积分,具体推导过程还需要经过傅里叶级数，不在此赘述）
$F(f)=\int^{\infty}_{-\infty}x(t)\mathrm{e}^{-\mathrm{i}(2\pi ft)}\mathrm{d}t \tag{10}$
$TC_k$ 的极限为 $F (f)$ ，根据式（5）得
$TC_k=\frac{T}{2}(A_k-\mathrm{i}B_k) \tag{11}$
根据式（8）可以得到，傅里叶谱振幅的最终结果要乘以 $T /2$

总结

python中scipy库的fft结果为 $y_k$ ，我们需要的傅里叶谱振幅为 $F_k$ ，根据式（5）（6）（7）（8）（11），得到二者的关系
$F_k = \frac{T}{2}X_k=\frac{T}{2}\cdot2|C_k|=\frac{T}{N}|y_k|=\Delta t\cdot|y_k|$
需要注意： $F_0 = \Delta t/2\cdot|y_0|$ ， $F_{N/2} = \Delta t/2\cdot|y_{N/2}|$