{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This R environment comes with many helpful analytics packages installed\n# It is defined by the kaggle/rstats Docker image: https://github.com/kaggle/docker-rstats\n# For example, here's a helpful package to load\n\nlibrary(tidyverse) # metapackage of all tidyverse packages\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nlist.files(path = \"../input\")\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2021-07-28T11:43:47.733960Z","iopub.execute_input":"2021-07-28T11:43:47.736132Z","iopub.status.idle":"2021-07-28T11:43:49.009425Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"install.packages(\"hdf5r\")","metadata":{"execution":{"iopub.status.busy":"2021-07-28T11:44:32.803986Z","iopub.execute_input":"2021-07-28T11:44:32.805514Z","iopub.status.idle":"2021-07-28T11:45:14.742713Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"library(RcppCNPy)\nlibrary(hdf5r)\nlibrary(signal)\nlibrary(tuneR)\nlibrary(seewave)\nlibrary(reticulate)\n","metadata":{"execution":{"iopub.status.busy":"2021-07-28T11:53:11.249819Z","iopub.execute_input":"2021-07-28T11:53:11.251569Z","iopub.status.idle":"2021-07-28T11:53:13.359785Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Actual GW sample data (from https://www.gw-openscience.org)\nfs = 4096  #sampling rate\n\nh1_file = H5File$new(\"../input/g2net-ligo-event-tutorial/LOSC_Event_tutorial/H-H1_LOSC_4_V2-1126259446-32.hdf5\", mode=\"r\")\nh1_time = c(1126259446.0, 1126259477.9997559)\nh1_tevent = 1126259462.44\n\nh1_strain_all = h1_file[[\"strain\"]][[\"Strain\"]]$read()\n\nh1_idx_start = floor((h1_tevent - 5 - h1_time[1])/(h1_time[2]-h1_time[1])*length(h1_strain_all))\nh1_idx_end = floor((h1_tevent + 5 - h1_time[1])/(h1_time[2]-h1_time[1])*length(h1_strain_all))\nh1_strain = h1_strain_all[h1_idx_start:h1_idx_end]\nh1_freq = seq(0, fs, by=fs/(length(h1_strain)-1))\n\nh1_file$close_all()\n","metadata":{"execution":{"iopub.status.busy":"2021-07-28T11:53:54.793877Z","iopub.execute_input":"2021-07-28T11:53:54.795290Z","iopub.status.idle":"2021-07-28T11:53:55.135587Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Whitening\nh1_fft = fft(h1_strain)\nplot(h1_freq, Re(h1_fft), type=\"l\", col=\"blue\", \n     xlab=\"Frequency\",\n     ylab=\"Amplitude\",\n     main=\"FFT of Actual sample data - H-H1_LOSC_4_V2-1126259446-32\")\n\nh1_asd_interp = approx(x=seq(0,length(h1_fft)-1), y=Mod(h1_fft), xout=h1_freq)\nnorm = 1/sqrt(1/(h1_freq[2]-h1_freq[1]))\n\n#h1_ifft = fft(h1_fft/abs(h1_fft)*norm, inverse = TRUE)/length(h1_fft)\nh1_ifft = fft(h1_fft/h1_asd_interp$y*norm, inverse = TRUE)/length(h1_fft)\n#h1_ifft = fft(h1_fft, inverse = TRUE)/length(h1_strain)\n#data1_ifft %>% Mod() %>% plot(log=\"xy\",  type=\"l\")\n\nplot(Re(h1_ifft), type=\"l\", col=\"blue\",\n     xlim=c(18000, 22000), ylim=c(-0.3, 0.3),\n     xlab=\"index\",\n     ylab=\"Strain\",\n     main=\"Whitened actual sample data - H-H1_LOSC_4_V2-1126259446-32\")\n","metadata":{"execution":{"iopub.status.busy":"2021-07-28T11:54:09.169462Z","iopub.execute_input":"2021-07-28T11:54:09.170988Z","iopub.status.idle":"2021-07-28T11:54:13.659371Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Filtering\nfn = fs/2  #Nyquist frequency\nfil_N = 4 #filter degree\nfc = c(43,800) #pass band frequency\nfc_norm = fc/fn  #Normalized frequency\n\nbutter_filt = butter(fil_N, fc_norm, type=\"pass\", plane=\"z\")\nfreqz(butter_filt, Fs=fs)  #Plot filter\n\nh1_filtered = filtfilt(butter_filt, h1_ifft)\nplot(h1_filtered, type=\"l\", xlim=c(3500,4000), ylim=c(-0.02, 0.02))\n","metadata":{"execution":{"iopub.status.busy":"2021-07-28T11:54:32.812978Z","iopub.execute_input":"2021-07-28T11:54:32.815184Z","iopub.status.idle":"2021-07-28T11:54:33.025683Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Spectrogram - all\n# We can find GW signal at the bottom center\nWave(left = h1_filtered, samp.rate = fs, bit = 16) %>%\n  spectro(osc = TRUE, \n          flim = c(0.03, fn/1000), \n          wl = fs/16, \n          collevels = seq(-110, 0, 0.1),\n#          flog = TRUE,\n          wn = \"blackman\",\n          palette = viridis::viridis,\n          ovlp = 15/16*100)\n\n# Spectrogram - close-up\n# Compliant with inspiral graviational wave form\nWave(left = h1_filtered, samp.rate = fs, bit = 16) %>%\n  spectro(osc = TRUE, \n          flim = c(0.03, 0.5), \n          wl = fs/16, \n          collevels = seq(-110, 0, 0.1),\n          tlim=c(4.5, 5.5),\n#          flog = TRUE,\n          wn = \"blackman\",\n          palette = viridis::viridis,\n          ovlp = 15/16*100)\n","metadata":{"execution":{"iopub.status.busy":"2021-07-28T11:54:41.340238Z","iopub.execute_input":"2021-07-28T11:54:41.342198Z","iopub.status.idle":"2021-07-28T11:56:24.040390Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Apply same approach to train data sample\nfs1 = 2048  #sampling rate\nfn1 = fs1/2  #Nyquist frequency\nfc1 = c(43,800) #pass band frequency\nfc1_norm = fc1/fn1  #Normalized frequency\nfil1_N = 4 #filter degree\n\n\nstrain = npyLoad(\"../input/g2net-gravitational-wave-detection/train/0/0/0/00000e74ad.npy\")[1,]\nfreq = seq(0, fs1, by=fs1/(length(strain)-1))\n","metadata":{"execution":{"iopub.status.busy":"2021-07-28T11:59:55.314230Z","iopub.execute_input":"2021-07-28T11:59:55.315997Z","iopub.status.idle":"2021-07-28T11:59:55.384114Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Whitening\ndata_fft = fft(strain)\nplot(freq, Re(data_fft), type=\"l\", col=\"blue\", \n     xlab=\"Frequency\",\n     ylab=\"Amplitude\",\n     main=\"FFT of sample data\")\n\nasd_interp = approx(x=seq(0,length(data_fft)-1), y=Mod(data_fft), xout=freq)\nnorm1 = 1/sqrt(1/(freq[2]-freq[1]))\n\ndata_ifft = fft(data_fft/asd_interp$y*norm1, inverse = TRUE)/length(data_fft)\n\n\n# Filtering\nbutter_filt1 = butter(fil1_N, fc1_norm, type=\"pass\", plane=\"z\")\ndata_filtered = filtfilt(butter_filt1, data_ifft)\n","metadata":{"execution":{"iopub.status.busy":"2021-07-28T12:00:19.587057Z","iopub.execute_input":"2021-07-28T12:00:19.588780Z","iopub.status.idle":"2021-07-28T12:00:19.699758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Spectrogram - all\n# Doesn't work as opposed to actual GW signal...\nWave(left = data_filtered, samp.rate = fs1, bit = 16) %>%\n  spectro(osc = TRUE, \n          flim = c(0.03, fn1/1000), \n          wl = fs1/16, \n          collevels = seq(-120, 0, 0.1),\n          #          flog = TRUE,\n          wn = \"blackman\",\n          palette = viridis::viridis,\n          ovlp = 15/16*100)\n\n","metadata":{"execution":{"iopub.status.busy":"2021-07-28T12:00:23.910533Z","iopub.execute_input":"2021-07-28T12:00:23.912002Z","iopub.status.idle":"2021-07-28T12:00:39.081915Z"},"trusted":true},"execution_count":null,"outputs":[]}]}