引力波的绘制

import numpy as np
import matplotlib.pyplot as plt  # 导入绘图使用的pyplot库
from scipy.io import wavfile  # 导入读取波形文件的scipy库

'''
可以使用wave.open(file, mode=None)来读.wav音频文件。
在打开文件的时候open(r'c:\')，加r和不加r是有区别的。
在字符串赋值的时候，前面加'r'可以防止字符串在时候的时候不被转义，原理是在转义字符前加'\'。
mode 可以为以下值：'rb'只读模式。'wb'只写模式。
它会返回一个元组，第一项为音频的采样率，第二项为音频数据的numpy数组。
wf_template.txt为引力波的理论模型，reftime为时间序列，ref_H1为信号数据。
genfromtxt主要执行两个循环运算。第一个循环将文件的每一行转换成字符串序列。
第二个循环将每个字符串序列转换为相应的数据类型。
transpose()转置。

'''
rate_h, hstrain = wavfile.read(r"H1_Strain.wav", "rb")
rate_l, lstrain = wavfile.read(r"L1_Strain.wav", "rb")
# reftime, ref_H1 = np.genfromtxt('GW150914_4_NR_waveform_template.txt').transpose()
reftime, ref_H1 = np.genfromtxt('wf_template.txt').transpose()  # 使用python123.io下载文件

'''
采样率的倒数为采样间隔。
plt.figure(figsize=(a,b))表示figure的大小为长、宽（单位为inch）
hstrain.shape[0]读取矩阵第一维度的长度，即数据点的个数，乘以采样间隔为坐标轴总长度。
要创造以原点为中心的序列，因此以-htime_len / 2为起点，htime_len / 2为终点，htime_interval为时间间隔。

'''
htime_interval = 1 / rate_h
ltime_interval = 1 / rate_l
fig = plt.figure(figsize=(12, 6))

# 丢失信号起始点
htime_len = hstrain.shape[0] / rate_h
htime = np.arange(-htime_len / 2, htime_len / 2, htime_interval)
plth = fig.add_subplot(221)
plth.plot(htime, hstrain, 'y')
plth.set_xlabel('Time (seconds)')
plth.set_ylabel('H1 Strain')
plth.set_title('H1 Strain')

ltime_len = lstrain.shape[0] / rate_l
ltime = np.arange(-ltime_len / 2, ltime_len / 2, ltime_interval)
pltl = fig.add_subplot(222)
pltl.plot(ltime, lstrain, 'g')
pltl.set_xlabel('Time (seconds)')
pltl.set_ylabel('L1 Strain')
pltl.set_title('L1 Strain')

pltref = fig.add_subplot(212)
pltref.plot(reftime, ref_H1)
pltref.set_xlabel('Time (seconds)')
pltref.set_ylabel('Template Strain')
pltref.set_title('Template')
fig.tight_layout()  # 自动调整图像外部边缘

plt.savefig("Gravitational_Waves_Original.png")
plt.show()
plt.close(fig)