地震子波模拟及其频谱分析实验

目的：
通过编写和运行相关程序，加深对地震子波类型和特征的理解。
内容及要求：
1、编写程序模拟及显示雷克子波，在同一张图上绘制主频为20Hz、35Hz、50Hz的零相位雷克子波波形。
2、分别绘制上述不同频率子波的振幅谱，对比振幅谱最大值对应频率与子波主频关系。

实验原理

时间域雷克子波函数为： ，其中为主频。对雷克子波频谱分析：对不同频率的雷克子波进行傅里叶变换得，所以不同频率的雷克子波振幅谱极大值对应频率即为主频。
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
代码如下
.
.
.
.
.
.
.
1.模拟及显示雷克子波Python代码

/**
 仅供学习，如有纰漏，望请指正
*/
1.import numpy as np
2.import matplotlib.pyplot as plt
3.import scipy
4.import math
5.
6.f1 = 20
7.f2 = 35
8.f3 = 50
9.
10.t = np.linspace(-0.055, 0.055, 500)
11.
12.rick1 = (1-2*(math.pi*f1*t)**2)*math.e**(-1*(math.pi*f1*t)**2)
13.
14.
15.rick2 = (1-2*(math.pi*f2*t)**2)*math.e**(-1*(math.pi*f2*t)**2)
16.
17.rick3 = (1-2*(math.pi*f3*t)**2)*math.e**(-1*(math.pi*f3*t)**2)
18.
19.plt.figure(figsize=(10,6))
20.fig,ax = plt.subplots()
21.ax.plot(t, rick1, '--', label="20Hz")
22.ax.plot(t, rick2, label="35Hz")
23.ax.plot(t, rick3, ':', label="50Hz")
24.ax.set_xlabel("Time")
25.ax.set_ylabel("Amp")
26.ax.set_title('Zero phase Rayleigh wavelet waveform')
27.ax.legend()
28.plt.show()

2.上述不同频率子波的振幅谱Python代码

/**
 仅供学习，如有纰漏，望请指正 --by wxw
*/
1.import numpy as np
2.import matplotlib.pyplot as plt
3.import scipy
4.import math
5.
6.f1 = 20
7.f2 = 35
8.f3 = 50
9.
10.f = np.linspace(0, 400, 500)
11.
12.R1 =(2*f**2)/(f1**2*f1*np.sqrt(math.pi))*math.e**(-1*(f/f1)**2)
13.R2 =(2*f**2)/(f2**2*f2*np.sqrt(math.pi))*math.e**(-1*(f/f2)**2)
14.R3 =(2*f**2)/(f3**2*f3*np.sqrt(math.pi))*math.e**(-1*(f/f3)**2)
15.
16.plt.figure(figsize=(10,6))
17.fig,ax = plt.subplots()
18.ax.plot(f, R1, '--', label="20Hz")
19.ax.plot(f, R2, label="35Hz")
20.ax.plot(f, R3, ':', label="50Hz")
21.ax.set_xlabel("Frequency")
22.ax.set_ylabel("Amp")
23.ax.set_title('Amplitude spectrum')
24.ax.legend()
25.plt.show()