#!/usr/bin/env python # coding: utf-8 # In[1]: import matplotlib import matplotlib.pyplot as plt import numpy as np from lib.PSD.LISAcalPSD import LISAcalPSD1, LISAcalPSD2 from lib.PSD.DECIGOcalPSD import B_DECIGOcalPSD, DECIGOcalPSD from pycbc.types.frequencyseries import FrequencySeries import pycbc.psd from scipy.interpolate import interp1d # interplate from scipy import signal # # 参数设置 # In[2]: sampFreq1 = 0.25 #采样频率(Sampling frequency)，单位时间样本点个数，应大于 2f（即Nyquist频率) duration1 = 2**16 #信号持续时间(duration of signal) n1 = int(duration1 * sampFreq1)#采样点数(Sampling Number), 有时也称为信号长度(Length of Signal) #为2的幂时，快速傅里叶变化效率最高 #n = duration * sampFreqint = (duration / sampIntrvl) sampIntrvl1 = 1.0 / sampFreq1 #采样周期(Sampling period)，隔多少时间取样一次，或步长 freqIntrvl1 = sampFreq1 / n1 #傅里叶变换 频率分辨率(Frequency Interval) # freqIntrvl = 1 / duration = 1 / (n * sampIntrvl) # = sampFreq / n f_max1 = sampFreq1/2.0 #信号模式的最大频率 print("采样频率为%3.2fHz，信号持续时间%ds, 时域信号采样%d 个点"%(sampFreq1, duration1, n1)) print("信号中可分析最大频率为%fHz"%f_max1) print("\n采样周期，即时域分辨率为%fs"%(sampIntrvl1)) print("信号频域的频率间隔，即频域分辨率为%fHz"%freqIntrvl1); # # 构建用于pyCBC的PSD频率序列 # ## 利用psd模拟函数生成 # In[3]: freqVec1 = np.linspace(0, freqIntrvl1 * (n1-1), n1) low_freq_cutoff1 = 1e-4 #构建频率序列(必须等频率间隔)，这里我们开始于0Hz(必须开始于0，否则导入pyCBC的结果不对), 结束于1Hz，频率间隔为上面的freqIntrvl #LISA PSD 模拟函数 psd1 = LISAcalPSD1(freqVec1) ''' 1. Babak, S., Fang, H., Gair, J. R., Glampedakis, K. & Hughes, S. A. Kludge’ gravitational waveforms for a test-body orbiting a Kerr black hole. Phys.Rev. D 75, 024005 (2007). ''' psd2 = LISAcalPSD2(freqVec1) ''' 1. Sathyaprakash, B. S. & Schutz, B. F. Physics, Astrophysics and Cosmology with Gravitational Waves. Living Rev Relativ 12, 122004 (2009). ''' #构建用于pyCBC的FrequencySeries psd1 = pycbc.psd.read.from_numpy_arrays(freqVec1 , psd1, n1, freqIntrvl1, low_freq_cutoff1) psd2 = pycbc.psd.read.from_numpy_arrays(freqVec1 , psd2, n1, freqIntrvl1, low_freq_cutoff1) # In[4]: # DECIGO sampFreq3 = 2**8 #采样频率(Sampling frequency)，单位时间样本点个数，应大于 2f（即Nyquist频率) duration3 = 2**14 #信号持续时间(duration of signal) n3 = int(duration3 * sampFreq3)#采样点数(Sampling Number), 有时也称为信号长度(Length of Signal) #为2的幂时，快速傅里叶变化效率最高 #n = duration * sampFreqint = (duration / sampIntrvl) sampIntrvl3 = 1.0 / sampFreq3 #采样周期(Sampling period)，隔多少时间取样一次，或步长 freqIntrvl3 = sampFreq3 / n3 #傅里叶变换 频率分辨率(Frequency Interval) # freqIntrvl = 1 / duration = 1 / (n * sampIntrvl) # = sampFreq / n f_max3 = sampFreq3/2.0 #信号模式的最大频率 print("采样频率为%3.2fHz，信号持续时间%ds, 时域信号采样%d 个点"%(sampFreq3, duration3, n3)) print("信号中可分析最大频率为%fHz"%f_max3) print("\n采样周期，即时域分辨率为%fs"%(sampIntrvl3)) print("信号频域的频率间隔，即频域分辨率为%fHz"%freqIntrvl3); freqVec3 = np.linspace(0, freqIntrvl3 * (n3-1), n3) low_freq_cutoff3 = 1e-2 #构建频率序列(必须等频率间隔)，这里我们开始于0Hz(必须开始于0，否则导入pyCBC的结果不对), 结束于1Hz，频率间隔为上面的freqIntrvl #LISA PSD 模拟函数 psd5 = B_DECIGOcalPSD(freqVec3) psd6 = DECIGOcalPSD(freqVec3) ''' 1. Multiband Gravitational-Wave Astronomy:Observing binary inspirals with a decihertz detector, B-DECIGO https://arxiv.org/pdf/1802.06977.pdf ''' #构建用于pyCBC的FrequencySeries psd5 = pycbc.psd.read.from_numpy_arrays(freqVec3 , psd5, n3, freqIntrvl3, low_freq_cutoff3) psd6 = pycbc.psd.read.from_numpy_arrays(freqVec3 , psd6, n3, freqIntrvl3, low_freq_cutoff3) # ## 从文件中导入 # In[5]: # 导入的psd文件要求第一列为频率(必须从0开始，不是等间隔需要进行插值)，第二列是psd，两列用空格分隔 ''' from WWW Sensitivity Curve Generator located at: http://www.srl.caltech.edu/~shane/sensitivity/MakeCurve.html ''' #导入 psdfile = np.fromfile("lib/PSD/LISAshotPSD.txt", dtype=float, count=-1, sep='\n') psdlen = len(psdfile) freqVecPre = psdfile[0:psdlen:2] psdPre = psdfile[1:psdlen:2] #插值 LISAcalPSD3 = interp1d(freqVecPre, psdPre, kind='cubic', fill_value="extrapolate") #构建用于pyCBC的FrequencySeries psd3 = LISAcalPSD3(freqVec1) psd3 = FrequencySeries(psd3, delta_f=freqIntrvl1, epoch='', dtype=None, copy=True) psd3 = pycbc.psd.read.from_numpy_arrays(freqVec1, psd3, n1, freqIntrvl1, low_freq_cutoff1) # ## pycbc内置 # In[6]: # 列出lalsuite内置的解析psd (没发现有LISA的，下面以LIGO的作为示例) print(pycbc.psd.get_lalsim_psd_list()) sampFreq2 = 2048.0 #采样频率(Sampling frequency)，单位时间样本点个数，应大于 2f（即Nyquist频率) duration2 = 4.0 #信号持续时间(duration of signal) n2 = int(duration2 * sampFreq2)#采样点数(Sampling Number), 有时也称为信号长度(Length of Signal) sampIntrvl2 = 1.0 / sampFreq2 #采样周期(Sampling period)，隔多少时间取样一次，或步长 freqIntrvl2 = sampFreq2 / n2 #傅里叶变换 频率分辨率(Frequency Interval) # freqIntrvl = 1 / duration = 1 / (n * sampIntrvl) # = sampFreq / n low_frequency_cutoff2 = 10.0 #低于此频率的psd将被设置为0 #示例，psd参见， https://dcc.ligo.org/LIGO-T1800044/public psd4 = pycbc.psd.from_string('aLIGOaLIGODesignSensitivityT1800044', n2, freqIntrvl2, low_frequency_cutoff2) # # 绘图 # In[7]: #绘制 频率 - sqrt(PSD) 图 #绘图字体参数设置 matplotlib.rcParams.update({'font.size': 15, 'font.family': 'STIXGeneral', 'mathtext.fontset': 'stix'}) plt.figure(figsize=(8,6)) plt.loglog(psd1.sample_frequencies[psd1>0], np.sqrt(psd1[psd1>0]), label = 'LISA_1') plt.loglog(psd2.sample_frequencies[psd2>0], np.sqrt(psd2[psd2>0]), label = 'LISA_2') plt.loglog(psd3.sample_frequencies[psd3>0], np.sqrt(psd3[psd3>0]), label = 'LISA_3') plt.loglog(psd4.sample_frequencies[psd4>0], np.sqrt(psd4[psd4>0]), label = 'LIGO') plt.loglog(psd5.sample_frequencies[psd5>0], np.sqrt(psd5[psd5>0]), label = 'B_DECIGO') plt.loglog(psd6.sample_frequencies[psd6>0], np.sqrt(psd6[psd6>0]), label = 'DECIGO') plt.xlim(1e-4, 2048) plt.xlabel("Frequency / Hz") plt.ylabel("$\sqrt{S_{n}(f) \ / \ (Hz^{-1})}$") plt.legend(loc = 'upper right') plt.grid(linestyle = "dotted", color = "#d3d3d3" , which="both") plt.tight_layout() plt.show()