import numpy as np import matplotlib.pyplot as plt import readligo as rl from scipy import signal from scipy.signal import butter, filtfilt, iirdesign, zpk2tf, freqz fn_L1 = "L-L1_LOSC_4_V2-1126259446-32.hdf5" strain_L1, time_L1, chan_dict_L1 = rl.loaddata(fn_L1, 'L1') # both H1 and L1 will have the same time vector, so: time = time_L1 # the time sample interval (uniformly sampled!) dt = time[1] - time[0] # Let's look at the data and print out some stuff: print( 'strain_L1: len, min, mean, max = ', \ len(strain_L1), strain_L1.min(),strain_L1.mean(),strain_L1.max()) #What's in chan_dict? (See also https://losc.ligo.org/tutorials/) bits = chan_dict_L1['DATA'] print("For L1, {0} out of {1} seconds contain usable DATA".format(bits.sum(), len(bits))) tevent = 1126259462.44 deltat = 5 indxt = np.where((time >= tevent-deltat) & (time < tevent+deltat)) plt.figure() plt.plot(time[indxt]-tevent,strain_L1[indxt],'g',label='L1 strain') plt.xlabel('time (s) since '+str(tevent)) plt.ylabel('strain') plt.legend(loc='lower right') plt.title('Advanced LIGO strain data') plt.savefig("strain_test.png") fs = 4096 from scipy.interpolate import interp1d import matplotlib.mlab as mlab make_plots = 1 make_psds = 1 if make_psds: # number of sample for the fast fourier transform: NFFT = 4*fs Pxx_L1, freqs = mlab.psd(strain_L1, Fs = fs, NFFT = NFFT) # We will use interpolations of the ASDs computed above for whitening: psd_L1 = interp1d(freqs, Pxx_L1) # Here is an approximate, smoothed PSD for H1 during O1, with no lines. We'll use it later. Pxx = (1.e-22*(18./(0.1+freqs))**2)**2+0.7e-23**2+((freqs/2000.)*4.e-23)**2 psd_smooth = interp1d(freqs, Pxx) if make_plots: # plot the ASDs, with the template overlaid: f_min = 20. f_max = 2000. plt.figure(figsize=(10,8)) plt.loglog(freqs, np.sqrt(Pxx_L1),'g',label='L1 strain') plt.loglog(freqs, np.sqrt(Pxx),'k',label='H1 strain, O1 smooth model') plt.axis([f_min, f_max, 1e-24, 1e-19]) plt.grid('on') plt.ylabel('ASD (strain/rtHz)') plt.xlabel('Freq (Hz)') plt.legend(loc='upper center') plt.title('Advanced LIGO strain data') plt.savefig("asd_test.png") """ from scipy.signal import butter, filtfilt fband=[43.0, 300.0] fs=4096 # We need to suppress the high frequency noise (no signal!) with some bandpassing: bb, ab = butter(4, [fband[0]*2./fs, fband[1]*2./fs], btype='band') strain_L1_bp = filtfilt(bb, ab, strain_L1) """ # generate linear time-domain filter coefficients, common to both H1 and L1. # First, define some functions: # This function will generate digital filter coefficients for bandstops (notches). # Understanding it requires some signal processing expertise, which we won't get into here. def iir_bandstops(fstops, fs, order=4): """ellip notch filter fstops is a list of entries of the form [frequency (Hz), df, df2] where df is the pass width and df2 is the stop width (narrower than the pass width). Use caution if passing more than one freq at a time, because the filter response might behave in ways you don't expect. """ nyq = 0.5 * fs # Zeros zd, poles pd, and gain kd for the digital filter zd = np.array([]) pd = np.array([]) kd = 1 # Notches for fstopData in fstops: fstop = fstopData[0] df = fstopData[1] df2 = fstopData[2] low = (fstop - df) / nyq high = (fstop + df) / nyq low2 = (fstop - df2) / nyq high2 = (fstop + df2) / nyq z, p, k = iirdesign([low,high], [low2,high2], gpass=1, gstop=6, ftype='ellip', output='zpk') zd = np.append(zd,z) pd = np.append(pd,p) # Set gain to one at 100 Hz...better not notch there bPrelim,aPrelim = zpk2tf(zd, pd, 1) outFreq, outg0 = freqz(bPrelim, aPrelim, 100/nyq) # Return the numerator and denominator of the digital filter b,a = zpk2tf(zd,pd,k) return b, a def get_filter_coefs(fs): # assemble the filter b,a coefficients: coefs = [] # bandpass filter parameters lowcut=43 highcut=300 order = 4 # bandpass filter coefficients nyq = 0.5*fs low = lowcut / nyq high = highcut / nyq bb, ab = butter(order, [low, high], btype='band') coefs.append((bb,ab)) # Frequencies of notches at known instrumental spectral line frequencies. # You can see these lines in the ASD above, so it is straightforward to make this list. notchesAbsolute = np.array( [14.0,34.70, 35.30, 35.90, 36.70, 37.30, 40.95, 60.00, 120.00, 179.99, 304.99, 331.49, 510.02, 1009.99]) # notch filter coefficients: for notchf in notchesAbsolute: bn, an = iir_bandstops(np.array([[notchf,1,0.1]]), fs, order=4) coefs.append((bn,an)) # Manually do a wider notch filter around 510 Hz etc. bn, an = iir_bandstops(np.array([[510,200,20]]), fs, order=4) coefs.append((bn, an)) # also notch out the forest of lines around 331.5 Hz bn, an = iir_bandstops(np.array([[331.5,10,1]]), fs, order=4) coefs.append((bn, an)) return coefs # and then define the filter function: def filter_data(data_in,coefs): data = data_in.copy() for coef in coefs: b,a = coef # filtfilt applies a linear filter twice, once forward and once backwards. # The combined filter has linear phase. data = filtfilt(b, a, data) return data fs = 4096. coefs = get_filter_coefs(fs) strain_L1_filt = filter_data(strain_L1, coefs) data = strain_L1_filt.copy() datsize = len(strain_L1) datafreq = np.fft.fftfreq(datsize)*fs dwindow = signal.tukey(datsize, alpha=1./8) data_fft = np.fft.fft(data*dwindow) / fs plt.figure(figsize=(15,15)) plt.subplot(4,1,1) plt.plot(datafreq,np.absolute(data_fft)) plt.axis([0,1200,0,1e-21]) data=strain_L1.copy() data_fft=np.fft.fft(data*dwindow)/fs ## test 2 Nc = 65536 filt_resp = np.ones(Nc) for coef in coefs: b,a = coef w,r = signal.freqz(b,a,worN=Nc) filt_resp = filt_resp*np.abs(r) freqf = (fs*.5 / np.pi) * w # We "double pass" the filtering using filtfilt, so we square the filter response filt_resp = filt_resp**2 print len(datafreq[:65536]), len(datafreq) plt.subplot(4,1,2) plt.plot(datafreq[:65536],np.absolute(data_fft)[:65536]) plt.axis([0,1200,0,1e-21]) plt.subplot(4,1,3) plt.plot(freqf,filt_resp) plt.axis([0,1200,0,1]) plt.subplot(4,1,4) plt.plot(datafreq[:65536],np.absolute(data_fft[65536:])*filt_resp,color="r") #plt.plot(freqf,filt_resp,color="g") plt.axis([20,1200,0,1e-21]) #plt.axis([20,1200,0,1]) plt.savefig("ft_test2.png")