{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"

🚀 G2Net Getting Started 🚀

\n\n***","metadata":{}},{"cell_type":"markdown","source":"
\n

\n Hopefully this notebook will give you a basic understanding of the task and data involved in this competition. Please give an upvote if you find it useful 👍\n

\n
\n \n
\"Space\"","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:40:27.984824Z","iopub.execute_input":"2022-10-26T21:40:27.985434Z","iopub.status.idle":"2022-10-26T21:40:27.991296Z","shell.execute_reply.started":"2022-10-26T21:40:27.985381Z","shell.execute_reply":"2022-10-26T21:40:27.990178Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import seaborn as sns","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:40:27.993471Z","iopub.execute_input":"2022-10-26T21:40:27.99385Z","iopub.status.idle":"2022-10-26T21:40:28.002319Z","shell.execute_reply.started":"2022-10-26T21:40:27.993804Z","shell.execute_reply":"2022-10-26T21:40:28.001316Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Contents:\n[Load in the data ⏳](#first-bullet)\n\n[What is HDF5 🤔](#second-bullet)\n \n[Plotting Spectograms 📊](#third-bullet) \n\n[What do the real and imaginary parts of a Short-time Fourier Transforms mean? 👻](#fourth-bullet) \n\n[Timestamp analysis ⏱](#fith-bullet)\n\n[Generating more data ⚙️](#sixth-bullet) \n\n[Baseline 📈](#seventh-bullry)","metadata":{"execution":{"iopub.status.busy":"2022-10-06T22:32:53.944249Z","iopub.execute_input":"2022-10-06T22:32:53.944766Z","iopub.status.idle":"2022-10-06T22:32:53.95639Z","shell.execute_reply.started":"2022-10-06T22:32:53.944726Z","shell.execute_reply":"2022-10-06T22:32:53.954342Z"}}},{"cell_type":"markdown","source":"### Task overview\n\n The goal of this competition is to find continuous gravitational-wave signals. You will develop a model sensitive enough to detect weak yet long-lasting signals emitted by rapidly-spinning neutron stars within noisy data. \n\n Each sample is comprised of a set of Short-time Fourier Transforms (SFTs) and corresponding time stamps for each interferometer. The SFTs are not always contiguous in time, since the interferometers are not continuously online.\n\n Each data sample contains either real or simulated noise and possibly a simulated continuous gravitational-wave signal (CW) . The task is to identify when a signal is present in the data (target=1).\n","metadata":{}},{"cell_type":"markdown","source":"\n### Continuous gravitational waves\nContinuous gravitational waves are produced by systems that have a __fairly constant and well-defined frequency__. Examples of these are binary sytems of stars or __black holes orbiting each other__ (long before merger), or a __single star swiftly rotating about its axis__ with a large mountain or other irregularity on it. These sources are expected to produce comparatively weak gravitational waves since they evolve over longer periods of time and are usually less catastrophic than sources producing inspiral or burst gravitational waves. The sound these gravitational waves would produce is a continuous tone since their frequency is nearly constant.\n> https://www.ligo.org/science/GW-Continuous.php\n\n\n\n\n\n \n\n\n\n \n\n\n\n \n \n
The Crab Nebula (catalogue designations M1, NGC 1952, Taurus A) is a supernova remnant and pulsar wind nebula in the constellation of Taurus. The common name comes from William Parsons, 3rd Earl of Rosse, who observed the object in 1842 using a 36-inch (91 cm) telescope and produced a drawing that looked somewhat like a crab. The nebula was discovered by English astronomer John Bevis in 1731, and it corresponds with a bright supernova recorded by Chinese astronomers in 1054. The nebula was the first astronomical object identified that corresponds with a historical supernova explosion.
First searches for gravitational waves from r-modes of the Crab pulsar... We did not find any evidence of gravitational waves. Our best 90% confidence upper limits on gravitational wave intrinsic strain were 1.5e-25 for the first run, 1.3e-25 for the first stretch of the second run, and 1.1e-25 for the second stretch of the second run... https://arxiv.org/abs/2101.00714
While LIGO and similar experiments detect objects in the final seconds before they merge, pulsar timing arrays are sensitive to gravitational wave signals from supermassive black holes that are spiraling toward each other and will not combine for millions of years. That's because galaxies merge hundreds of millions of years before the central black holes they host combine to make one giant supermassive black hole....\n \"A difference between when the pulsar signals should arrive, and when they do arrive, can signal a gravitational wave,\" Mingarelli said. \"And since the pulsars we study are about 3,000 light-years away, they act as a galactic-scale gravitational-wave detector.\" https://www.nasa.gov/feature/jpl/listening-for-gravitational-waves-using-pulsars
\n\n> When stars at least eight times as massive as our Sun run out of nuclear fuel, their lives end in supernovae. These explosions bring forth new cosmic objects: neutron stars. Because of their strong magnetic fields and fast rotation, neutron stars emit beamed radio waves and gamma rays like a cosmic lighthouse. And just like a sailor will see flashes of light from the lighthouse, terrestrials will see the neutron star as a pulsating radio or gamma-ray source: a pulsar.\n> Cosmic lighthouses born in supernova explosions\n> \n> Tiny bumps on neutron stars\n> Because of their strong gravity, neutron stars are expected to be almost perfectly spherical – every mountain would be flattened because regular matter cannot withstand the strong gravitational force at the surface. Perfectly round spheres isolated in space will not emit gravitational waves, no matter how fast they spin. If, however, a spinning neutron star is non-axisymmetrically deformed, it will emit gravitational waves.\n> \n> Those deformations can come into existence because the crystalline crust is anything but regular matter. It is 20 times as hard as steel and could thus support small “bumps” on the neutron star surface that are a few centimeters high. They could form in the turbulent aftermath of the supernova or when matter from a stellar companion piles up on the neutron star’s surface. The magnitude of neutron star deformations is usually measured as an “ellipticity”, which describes how much the neutron star deviates from a perfect sphere. In case of those bumps, the ellipticity can be as large as 10-6. Other deformations could result from the neutron star’s magnetic field – it is one billion times as strong as the Earth’s. Simulations of these magnetic deformations show that they can induce smaller ellipticities as large as 10-8.\n> \n> Independently of what causes the deformations: for as long they are there, a spinning neutron star will emit gravitational waves at __twice its rotation frequency__. This is similar to binary systems, where gravitational waves are emitted at twice the orbital frequency. About a fifth of the known neutron stars (which have been observed as pulsars) spin at frequencies of more than __5 Hertz (and up to 700 Hertz)__. Their expected gravitational-wave emission will therefore be in the frequency band observed by ground-based detectors such as LIGO.\n\n","metadata":{}},{"cell_type":"markdown","source":"### **Create clean GW signal (chirp)**\nmost algoritms study very short gravitational-wave signals, produced in the final seconds of black hole and neutron star coalescences. __This is not a continous \"ON\" gravitational wave.__\n\n> Strong and transient versus faint and continuous\n> So far, all gravitational-wave events observed by the LIGO and Virgo detectors came from merging compact binary systems. They observed the inspiral and merger of several dozens of black holes and/or neutron stars in binary systems. The gravitational waves emitted from these compact binary mergers are transient: They are in the LIGO/Virgo frequency band for fractions of a second to up to more than a minute before they end.\nIn contrast, the gravitational waves emitted by rapidly spinning neutron stars are continuous and much fainter than the transient signals from compact binary systems – at least __by a factor of 10,000__. On the positive side, __they are continuous, and the signal is simple__. By analyzing long stretches of data, even very faint signals can be pulled out of noisy detector data. Once a signal has been identified, it can be repeatedly observed as detectors and analysis methods improve. But the search is computationally very demanding.\n\n> https://www.einstein-online.info/en/spotlight/continuousGW/\n\nReverse engineering: Create clean GW signals -GEIR DRANGE\nRiroriro: Simulating gravitational waves and evaluating their detectability in Python\nRiroriro is a Python package to simulate the gravitational waveforms of __binary mergers\nof black holes and/or neutron stars__, and calculate several properties of these mergers and\nwaveforms, specifically relating to their observability by gravitational wave detectors.","metadata":{}},{"cell_type":"code","source":"!pip install riroriro","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:40:28.003521Z","iopub.execute_input":"2022-10-26T21:40:28.004347Z","iopub.status.idle":"2022-10-26T21:41:06.363019Z","shell.execute_reply.started":"2022-10-26T21:40:28.004299Z","shell.execute_reply":"2022-10-26T21:41:06.361127Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import riroriro.inspiralfuns as ins\nimport riroriro.mergerfirstfuns as me1\nimport riroriro.matchingfuns as mat\nimport riroriro.mergersecondfuns as me2\nimport librosa\nimport librosa.display","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:41:06.366348Z","iopub.execute_input":"2022-10-26T21:41:06.366808Z","iopub.status.idle":"2022-10-26T21:41:07.071227Z","shell.execute_reply.started":"2022-10-26T21:41:06.36676Z","shell.execute_reply":"2022-10-26T21:41:07.069729Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import math\nimport os\nimport matplotlib.pyplot as plt\nfrom ipywidgets import interact, fixed, FloatSlider\nimport IPython.display as ipd\n%matplotlib inline","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:41:07.073855Z","iopub.execute_input":"2022-10-26T21:41:07.074789Z","iopub.status.idle":"2022-10-26T21:41:07.085636Z","shell.execute_reply.started":"2022-10-26T21:41:07.074716Z","shell.execute_reply":"2022-10-26T21:41:07.083888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Parameters:\n# logMC: system mass (0.0-2.0)\n# q: mass ratio (0.1-1.0)\n# D: distance (Mpc)\n# merger_type: 'BH'=binary black hole merger, 'NS'=binary neutron star merger\n# flow: low frequency (Hz) \ndef gen_gw(logMc=1.4, q=0.8, D=100.0, flow=10.0, merger_type='BH'):\n M, eta = ins.get_M_and_eta(logMc=logMc,q=q)\n start_x = ins.startx(M,flow)\n end_x = ins.endx(eta,merger_type)\n x, xtimes, dt = ins.PN_parameter_integration(start_x,end_x,M,eta)\n realtimes = ins.inspiral_time_conversion(xtimes,M)\n i_phase, omega, freq = ins.inspiral_phase_freq_integration(x,dt,M)\n r, rdot = ins.radius_calculation(x,M,eta)\n A1, A2 = ins.a1_a2_calculation(r,rdot,omega,D,M,eta)\n i_Aorth, i_Adiag = ins.inspiral_strain_polarisations(A1,A2,i_phase)\n i_amp = ins.inspiral_strain_amplitude(i_Aorth,i_Adiag)\n i_time = realtimes\n i_omega = omega\n sfin, wqnm = me1.quasi_normal_modes(eta)\n alpha, b, C, kappa = me1.gIRS_coefficients(eta,sfin)\n fhat, m_omega = me1.merger_freq_calculation(wqnm,b,C,kappa)\n fhatdot = me1.fhat_differentiation(fhat)\n m_time = me1.merger_time_conversion(M)\n min_switch_ind = mat.min_switch_ind_finder(i_time,i_omega,m_time,m_omega)\n final_i_index = mat.final_i_index_finder(min_switch_ind,i_omega,m_omega)\n time_offset = mat.time_offset_finder(min_switch_ind,final_i_index,i_time,m_time)\n i_m_time, i_m_omega = mat.time_frequency_stitching(min_switch_ind,final_i_index,time_offset,i_time,i_omega,m_time,m_omega)\n i_m_freq = mat.frequency_SI_units(i_m_omega,M)\n m_phase = me2.merger_phase_calculation(min_switch_ind,final_i_index,i_phase,m_omega)\n i_m_phase = me2.phase_stitching(final_i_index,i_phase,m_phase)\n m_amp = me2.merger_strain_amplitude(min_switch_ind,final_i_index,alpha,i_amp,m_omega,fhat,fhatdot)\n i_m_amp = me2.amplitude_stitching(final_i_index,i_amp,m_amp)\n m_Aorth, m_Adiag = me2.merger_polarisations(final_i_index,m_amp,m_phase,i_Aorth)\n i_m_Aorth, i_m_Adiag = me2.polarisation_stitching(final_i_index,i_Aorth,i_Adiag,m_Aorth,m_Adiag)\n return np.array(i_m_time), np.array(i_m_Aorth), np.array(i_m_Adiag), np.array(i_m_freq)","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:41:07.087373Z","iopub.execute_input":"2022-10-26T21:41:07.087757Z","iopub.status.idle":"2022-10-26T21:41:07.103614Z","shell.execute_reply.started":"2022-10-26T21:41:07.087722Z","shell.execute_reply":"2022-10-26T21:41:07.102565Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The function returns two waves that represent orthogonal/diagonal waves. The output timescale that is returned is non-linear, so to convert these signals into uniform sampled signals as in the dataset, we need to resample. The function below will resample the gravitational wave signals to 2048Hz. It is crude though, based on nearest sample, but good enough for studying spectrums. Interpolation would be more proper","metadata":{}},{"cell_type":"code","source":"SR = 2048 # target sample rate (Hz)\n# Parameters:\n# dt: time series\n# amp: amplitude signal\n# seg: output sequence length (seconds)\ndef resample(dt, amp, seg=2.0):\n end = dt[-1]\n start = end - seg\n d = np.zeros(int(SR*seg))\n for i in range((int(SR*seg))):\n t = start + i/SR\n d[i] = amp[np.where(dt == dt[np.abs(dt-t).argmin()])[0][0]]\n return d","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:41:07.105149Z","iopub.execute_input":"2022-10-26T21:41:07.105541Z","iopub.status.idle":"2022-10-26T21:41:07.12114Z","shell.execute_reply.started":"2022-10-26T21:41:07.105507Z","shell.execute_reply":"2022-10-26T21:41:07.119859Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Another helper function for plotting the last part of the GW signal (containing the chirp):","metadata":{}},{"cell_type":"code","source":"def plot_sig(dt, sig1, sig2=None, seg=2.0):\n end = dt[-1]\n start = end - seg\n plt.figure(1)\n plt.plot(dt, sig1)\n peak = np.max(np.abs(sig1))\n plt.axis([start,end,np.min(sig1)-peak/10,np.max(sig1)+peak/10])\n if sig2 is not None:\n plt.plot(dt, sig2)\n plt.xlabel('Time (s)')\n plt.ylabel('Strain amplitude')","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:41:07.122912Z","iopub.execute_input":"2022-10-26T21:41:07.123369Z","iopub.status.idle":"2022-10-26T21:41:07.132885Z","shell.execute_reply.started":"2022-10-26T21:41:07.12333Z","shell.execute_reply":"2022-10-26T21:41:07.131685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Test the signal generation¶\nNow, let's test the code by emulating the first GW detected (GW150914).","metadata":{}},{"cell_type":"code","source":"m_time, m_Aorth, m_Adiag, m_freq = gen_gw(logMc=1.4, q=0.2)","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:41:07.134252Z","iopub.execute_input":"2022-10-26T21:41:07.134635Z","iopub.status.idle":"2022-10-26T21:42:44.779055Z","shell.execute_reply.started":"2022-10-26T21:41:07.134601Z","shell.execute_reply":"2022-10-26T21:42:44.777858Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(16,8))\nplt.subplot(2, 1, 1)\nplot_sig(m_time, m_Aorth, m_Adiag, seg=2)\nplt.subplot(2, 1, 2)\nplot_sig(m_time, m_Aorth, m_Adiag, seg=.1)","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:42:44.785138Z","iopub.execute_input":"2022-10-26T21:42:44.785557Z","iopub.status.idle":"2022-10-26T21:42:45.682408Z","shell.execute_reply.started":"2022-10-26T21:42:44.785521Z","shell.execute_reply":"2022-10-26T21:42:45.681527Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Changing parameters: \n* logMC: system mass (0.0-2.0)\n* q: mass ratio (0.1-1.0)\n* D: distance (Mpc)\n* merger_type: 'BH'=binary black hole merger, 'NS'=binary neutron star merger\n* flow: low frequency (Hz) \nsystem mass (0.0-2.0)\nmass =2\nmass ratio =1\n","metadata":{}},{"cell_type":"code","source":"m_time, m_Aorth, m_Adiag, m_freq = gen_gw(logMc=1.4, q=1.0)\nfig = plt.figure(figsize=(16,8))\nplt.subplot(2, 1, 1)\nplot_sig(m_time, m_Aorth, m_Adiag, seg=2)\nplt.subplot(2, 1, 2)\nplot_sig(m_time, m_Aorth, m_Adiag, seg=.1)","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:42:45.683798Z","iopub.execute_input":"2022-10-26T21:42:45.684412Z","iopub.status.idle":"2022-10-26T21:44:38.833833Z","shell.execute_reply.started":"2022-10-26T21:42:45.684375Z","shell.execute_reply":"2022-10-26T21:44:38.832426Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"'NS'=binary neutron star merger\n\nFor binary neutron star and black hole–neutron star mergers, use 'NS'. For a 'NS' waveform, you should only simulate the inspiral portion (inspiralfuns) and not the merger/ringdown (as the mathematical equations used for the merger/ringdown are only valid for black holes and for neutron star mergers this phase would be beyond the detectable range of LIGO anyway),","metadata":{}},{"cell_type":"markdown","source":"\nIn practice, __the most important__ LIGO search targets for continuous gravitational waves __are neutron stars__ (ultra-compact stellar remnants) in our own Galaxy. Once we reach sufficient detector sensitivity to find such a signal, we expect we can observe it continuously for many months or years.\nOver short times, a continuous gravitational-wave signal from a Galactic neutron star will look almost perfectly constant in both frequency and amplitude, as seen in Fig. 1 for a 0.1 second interval. However, over longer durations, the frequency of the signal will slowly change, for two reasons. The first reason is that, as the neutron star emits gravitational and electromagnetic waves, it loses energy which causes it to rotate more slowly. The second reason is that the detector here on Earth is moving with respect to the neutron star, which changes the frequency of the gravitational waves observed in the detector. The manner in which the signal frequency changes is illustrated in Fig. 2.\n> https://www.ligo.org/science/GW-Continuous.php\n\n\n\n\n\n
\n","metadata":{}},{"cell_type":"markdown","source":"The top panel shows the frequency, in blue, changing with a daily cycle due to the rotation of the Earth. The middle panel zooms out to show the frequency, in red, changing on a year's time scale due to the orbit of the Earth around the Sun. The bottom panel zooms out further to show how the frequency, in green, slowly decreases due to the rotation of the neutron star itself slowing down over many years.\n\nTracking all possible frequency changes is what makes the detection of continuous gravitational waves a computational challenge.\n","metadata":{}},{"cell_type":"markdown","source":"Gravitational waves are ripples in the curvature of space-time and manifest themselves as fluctuating tidal forces on masses in the path of the wave. The first gravitational-wave detectors were based on the effect of these forces on the fundamental resonant mode of aluminium bars at room temperature. Initial instruments were constructed by Joseph Weber [310, 311] and subsequently developed by others.\nhttps://link.springer.com/article/10.12942/lrr-2011-5\n\n\n\n\n\n
\n\nThese massive aluminium cylinders vibrated at a resonance frequency of 1660 hertz and were designed to be set in motion by gravitational waves predicted by Weber. Because these waves were supposed to be so weak, the cylinders had to be massive and the piezoelectric sensors had to be very sensitive, capable of detecting a change in the cylinders' lengths by about 10−16 meters.[1]\n\nHowever, the most promising design of gravitational-wave detectors, offering the possibility of very high sensitivities over a wide range of frequency, uses widely-separated test masses freely suspended as pendulums on Earth or in a drag-free craft in space; laser interferometry provides a means of sensing the motion of the masses produced as they interact with a gravitational wave.\n\nGravitational-wave detectors of long baseline have been built in a number of places around the world;\n\nDetection methods cover over 13 orders of magnitude in frequency (see Figure 1) equivalent to covering from radio waves to X-rays in the electromagnetic spectrum. This broadband coverage allows to probe a wide range of potential sources.\n\n\n\n\n \n
\n\n\n\n","metadata":{}},{"cell_type":"markdown","source":"\n\n","metadata":{}},{"cell_type":"markdown","source":"## Detection of Gravitational Waves\nGravitational waves are most simply thought of as ripples in the curvature of space-time, their effect being to change the separation of adjacent masses on Earth or in space; this tidal effect is the basis of all present detectors. Gravitational wave strengths are characterised by the gravitational-wave amplitude h, given by\n\n$$h = {{2\\Delta L} \\over L},$$\nwhere ΔL is the change in separation of two masses a distance L apart; for the strongest-allowed component of gravitational radiation, the value of h is proportional to the third time derivative of the quadrupole moment of the source of the radiation and inversely proportional to the distance to the source. The radiation field itself is quadrupole in nature and this shows up in the pattern of the interaction of the waves with matter.\n\nThe problem for the experimental physicist is that the predicted magnitudes of the amplitudes or strains in space in the vicinity of the Earth caused by gravitational waves even from the most violent astrophysical events are extremely small, of the order of 10−21 or lower.\n\nhttps://link.springer.com/article/10.12942/lrr-2011-5\n\nIndeed, current theoretical models on the event rate and strength of such events suggest that in order to detect a few events per year — from coalescing neutron-star binary systems, for example, an amplitude sensitivity close to 10−22 over timescales as short as a millisecond is required. If the Fourier transform of a likely signal is considered it is found that the energy of the signal is distributed over a frequency range or bandwidth, which is approximately equal to 1/timescale. For timescales of a millisecond the bandwidth is approximately 1000 Hz, and in this case the spectral density of the amplitude sensitivity is obtained by dividing 10−22 by the square root of 1000. Thus, detector noise levels must have an amplitude spectral density lower than ≃ 10−23 Hz−1/2 over the frequency range of the signal. Signal strengths at the Earth, integrated over appropriate time intervals, for a number of sources are shown in Figure 2.\n\n\n\n\n \n
\n\n\n\n\n \n
\n The sensitivity of an interferometric gravitational-wave detector is limited by noise from various sources. Taking this frequency-dependent noise floor into account, a design goal can be estimated for a particular detector design. For example, the design sensitivity for initial design of LIGO is show plotted alongside the achieved sensitivities of the three individual interferometers during the fifth science run. Such strain sensitivities are expected to allow a reasonable probability for detecting gravitational wave sources. However, in order to guarantee the observation of a full range of sources and to initiate gravitational-wave astronomy, a sensitivity or noise performance approximately ten times better in the mid-frequency range and several orders of magnitude better at 10 Hz, is desired. \n
\n\n\n\n\n \n
\n","metadata":{}},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"### Main Noise Sources\n#### Seismic noise\nSeismic noise at a reasonably quiet site on the Earth follows a spectrum in all three dimensions close to 10−7f−2 m/Hz1/2 (where here and elsewhere we measure f in Hz) and thus if the disturbance to each test mass must be less than 3 × 10−20 m/Hz1/2 at, for example, 30 Hz, then the reduction of seismic noise required at that frequency in the horizontal direction is greater than 10E9. \n#### Gravity gradient (Newtonian) noise\nGravity gradients, caused by direct gravitational coupling of mass density fluctuations to the suspended mirrors, were identified as a potential source of noise in ground-based gravitational-wave detectors in 1972\n#### Thermal noise\nThermal noise associated with the mirror masses and the last stage of their suspensions is the most significant noise source at the low frequency end of the operating range of initial long baseline gravitational wave detectors\n#### Quantum noise\nPhotoelectron shot noise\nFor gravitational-wave signals to be detected, the output of the interferometer must be held at one of a number of possible points on an interference fringe\n#### Radiation pressure noise\nAs the effective laser power in the arms is increased, another phenomenon becomes increasingly important arising from the effect on the test masses of fluctuations in the radiation pressure. One interpretation on the origin of this radiation pressure noise may be attributed to the statistical uncertainty in how the beamsplitter divides up the photons of laser light","metadata":{}},{"cell_type":"markdown","source":"### Gravitational-wave sensitivity curves\n> https://cplberry.com/2015/01/10/1408-0740/\n\n> When talking about gravitational-wave detectors, you normally use a sensitivity curve. This shows how sensitive it is at a given frequency: you plot a graph with the sensitivity curve on, and then plot the spectrum of the source you’re interested in on the same graph. If your source is above the sensitivity curve, you can detect it (yay), but if it lies below it, then you can’t pick it out from the noise (boo). Making a plot with lots of sensitivity curves on sounds simple: you look up the details for lots of detectors, draw them together and add a few sources.\n\n\n\n\n \n \n \n\n \n \n \n \n\n \n \n \n
Gravitational-wave sensitivity-curve plot using characteristic strain. The area between the detector’s curve and the top of the box for a source indicates how loud that signal would be.
Gravitational-wave sensitivity-curve plot using the square root of the power spectral density (the amplitude spectral density)
Gravitational-wave sensitivity-curve plot using the energy density that cosmologists love. The proper name of the plotted quantity is the critical energy density per logarithmic frequency interval multiplied by the reduced Hubble constant squared.)
\n\nGravitational-wave sensitivity-curve plot using characteristic strain. The area between the detector’s curve and the top of the box for a source indicates how loud that signal would be.","metadata":{}},{"cell_type":"markdown","source":"### A milestone in the race to find the first continuous gravitational wave signal\n\n> https://www.aei.mpg.de/581560/a-milestone-in-the-race-to-find-the-first-continuous-gravitational-wave-signal\n\n> All so far observed gravitational-wave signals have been very short, produced in the final seconds of black hole and neutron star coalescences. There exist other signals, that are continously “ON”. Such continuous signals are much weaker than the short-lived ones but this drawback can be partially made up by utilizing more data and at the cost of a much greater algorithmic complexity of the search. One promising source of continuous gravitational waves are slightly deformed (non-axisymmetric) neutron stars. Our new results explore deformations smaller than 10 parts in a billion. This corresponds to star deformations as large as several widths of human hair! We find no signal candidates but identify a few outliers. Albeit it is unlikely that these candidates come from signals that precisely follow the assumed model, and may well come from coherent disturbances in the data, they will be explored very carefully as new data becomes available – something one can only do with continuous waves. Between the two papers, we search for signals from anywhere in the sky and with frequencies between 20 Hz and 1700 Hz.\n\nWe present the results of an all-sky search for continuous gravitational wave signals with frequencies in the __500-1700 Hz range__ targeting neutron stars with ellipticity of 1e-8. The search is done on LIGO O2 data using the Falcon analysis pipeline. The results presented here double the sensitivity over any other result on the same data. The search is capable of detecting low ellipticity sources up to 170 pc. We establish strict upper limits which hold for worst-case signal parameters. We list outliers uncovered by the search, including several which we cannot associate with any known instrumental cause.\n\n\n\n\n \n \n
Upper limits on the ellipticity of a source at a certain distance (black). We also show the recent upper limits from the low ellipticity all-sky search of Dergachev & Papa (2020). The dashed line is the spin-down ellipticity for the highest spin-down rate probed by each search.
\n","metadata":{}},{"cell_type":"markdown","source":"### Short-time Fourier Transforms (SFTs)\n*a wikipedia brief\nShort-time Fourier transform (STFT), is a Fourier-related transform used to determine the sinusoidal frequency and phase content of local sections of a signal as it changes over time.\n\nIn practice, the procedure for computing STFTs is to divide a longer time signal into shorter segments of equal length and then compute the Fourier transform separately on each shorter segment. This reveals the Fourier spectrum on each shorter segment. One then usually plots the changing spectra as a function of time, known as a spectrogram or waterfall plot, such as commonly used in software defined radio (SDR) based spectrum displays. \nFull bandwidth displays covering the whole range of an SDR commonly use fast Fourier transforms (FFTs) with 2^24 points on desktop computers.\n\nUsing the following sample signal $x(t)$that is composed of a set of four sinusoidal waveforms joined together in sequence. Each waveform is only composed of one of four frequencies (10, 25, 50, 100 Hz). The definition of $x(t)$ is:\n$$x(t)=\\begin{cases}\n\\cos (2 \\pi 10 t) & 0\\,\\mathrm{s} \\le t < 5 \\,\\mathrm{s} \\\\\n\\cos (2 \\pi 25 t) & 5\\,\\mathrm{s} \\le t < 10\\,\\mathrm{s} \\\\\n\\cos (2 \\pi 50 t) & 10\\,\\mathrm{s} \\le t < 15\\,\\mathrm{s} \\\\\n\\cos (2 \\pi 100 t) & 15\\,\\mathrm{s} \\le t < 20\\,\\mathrm{s} \\\\\n\\end{cases}$$\n\nThen it is sampled at 400 Hz. The following spectrograms were produced:\n\n\n\n\n\n\n\n\n\n
\n\n\n \n \n \n
\n The 25 ms window allows us to identify a precise time at which the signals change but the precise frequencies are difficult to identify. At the other end of the scale, the 1000 ms window allows the frequencies to be precisely seen but the time between frequency changes is blurred.\n
\n ","metadata":{}},{"cell_type":"markdown","source":"The Short Time Fourier Transform (STFT) is a special flavor of a Fourier transform where you can see how your frequencies in your signal change through time. It works by slicing up your signal into many small segments and taking the fourier transform of each of these. The result is usually a waterfall plot which shows frequency against time. \n\n\n\n\n \n \n \n\n \n \n \n
s a spectrum of a whistle
Here is a picture of a whistle which increases in frequency over time.It goes from 880 Hz (A in the 5th octave) all the way up to 1760 Hz (A in the 6th octave).
","metadata":{}},{"cell_type":"markdown","source":"#### STFT Algorithm:\n\nSo, we understand what we’re trying to make – now we have to figure out how to make it. The data flow we have to achieve is pretty simple, as we only need to do the following steps:\n\nPick out a short segment of data from the overall signal\nMultiply that segment against a half-cosine function\nPad the end of the segment with zeros\nTake the Fourier transform of that segment and normalize it into positive and negative frequencies\nCombine the energy from the positive and negative frequencies together, and display the one-sided spectrum\nScale the resulting spectrum into dB for easier viewing\nClip the signal to remove noise past the noise floor which we don’t care about\n\n* Step 1 – Pick segment: We need to find our current segment to process from the overall data set. We use the concept of a ‘sliding window’ to help us visualize what’s happening. The data inside the window is the current segment to be processed.\n\n\n\n\n\n\n \n \n \n\n\n \n \n \n\n\n \n \n \n\n\n \n \n \n\n \n \n \n \n\n \n \n \n \n\n \n \n \n \n\n \n \n
The window’s length remains the same during the processing of the data, but the offset changes with each step of the algorithm. Usually when processing the STFT, the change in offset will be less than one window length, meaning that the last window and the current window overlap. If we define the window size, and the percentage of overlap, we know all the information we need about how the window moves throughout the processing.
Step 2– Multiply by half-cosine: This helps to alleviate a problem created by segmenting the data. When the data is cut into pieces, the edges make a sharp transition that didn’t exist before. Multiplying by a half cosine function helps to fade the signal in and out so that the transitions at the edges do not affect the Fourier transform of the data.\n\n.
There is an additional benefit to using a half cosine window. It can be shown that with an overlap factor of 50%, the sum of the windows is a constant value of 1. This means that multiplying our data by this particular function does not introduce differences in amplitude from the original signal. This is mostly important if the waveform has some type of units attached to it [pascals] and we want the frequencies in the Fourier transform to correctly represent the energy content of each frequency in those units.
Step 3 – Pad Data: We pad the end of the current segment with a number of zeros equal to the length of the window. Or in other words, the new segment will be twice as long as the original segment. The reason for this is explained later.\n\n.
Step 4 – Fourier Transform and Scale: This is where all the magic happens. Everything up to this point was just preparing the data, but this step is the beating heart of the algorithm. Here we finally take the Fourier transform of the data. I’m not going to explain the Fourier transform here except to say that it takes a signal and extracts the frequency content from that signal. There are a few other things that need to be remembered about the Fourier transform.\n\nThe number of samples of data in is equal to the number of frequency bins out\nThe maximum frequency that can be represented is the Nyquist frequency, or half the sampling frequency of the data.\nIn the the Discrete Fourier Transform (what we are using) the order of the frequency bins is 0 Hz (DC component), positive frequencies, Nyquist Frequency, negative frequencies\nThe data from the Fourier transform needs to be scaled by the number of samples in the transform to maintain equal energy (according to Parseval’s theorem)\n\n.
Step 5 – Autopower: This is just more post processing to make the results more palatable to look at. We find the autopower spectrum from the results of the Fourier transform. The autopower spectrum is a one-sided spectrum (only contains positive frequencies) as opposed to the two-sided spectrum returned by the Fourier transform – so it’s more like the actual real world which does not contain negative frequencies. The autopower spectrum is also scaled to represent the energy content that was in both the positive and negative frequencies of the Fourier transform. This is important for the ‘correctness’ of the energy content of each frequency, but also has the handy benefit of suppressing lower energy content and dramatizing higher energy content. Usually people only care about the higher energy content, so finding this transform typically makes the data much easier to look at.\n\n.
Step 6 – Convert to DB: This is another transform to make the data easier to look at. It turns out that peoples ears work on a logarithmic scale so that the ear can detect much finer changes in amplitude at low amplitudes than at high amplitudes. We transform the data into decibels (dB) which is a logarithmic scale, so that we can see the energy content of a signal more how our ears would detect it. Converting the data to decibels has the effect of stretching peaks downwards towards the average sound level, and bringing troughs upwards. This allows us to compare content at all amplitude levels. If we didn’t do this and we scaled the data so that the peaks were visible, it would be impossible to see any shape in the quieter values, vs if we looked at the quieter values, the peaks would explode off the top of the screen.\n\n.
Step 7 – Clip Data: This also makes the data easier to look at. We know from experience that everything below -40 dB is well below the noise floor and is probably just numerical error in the algorithm. Therefore, we can clip the data so that everything below -40 dB is set to -40 dB exactly. This gives us more color range to apply to significant portions of our data.\n.
\n\n","metadata":{}},{"cell_type":"markdown","source":"> The Fourier transform yields frequency information that is averaged over the entire time domain. However, the information on when these frequencies occur is hidden in the transform. This phenomenon is illustrated by the following example.\n\n> https://www.audiolabs-erlangen.de/resources/MIR/FMP/C2/C2_STFT-Basic.html","metadata":{}},{"cell_type":"code","source":"\n\nFs = 128\nduration = 10\nomega1 = 1\nomega2 = 5\nN = int(duration * Fs)\nt = np.arange(N) / Fs\nt1 = t[:N//2]\nt2 = t[N//2:]\n\nx1 = 1.0 * np.sin(2 * np.pi * omega1 * t1)\nx2 = 0.7 * np.sin(2 * np.pi * omega2 * t2)\nx = np.concatenate((x1, x2))\n\nplt.figure(figsize=(8, 2))\nplt.subplot(1, 2, 1)\nplt.plot(t, x, c='k')\nplt.xlim([min(t), max(t)])\nplt.xlabel('Time (seconds)')\n\nplt.subplot(1, 2, 2)\nX = np.abs(np.fft.fft(x)) / Fs\nfreq = np.fft.fftfreq(N, d=1/Fs)\nX = X[:N//2]\nfreq = freq[:N//2]\nplt.plot(freq, X, c='k')\nplt.xlim([0, 7])\nplt.ylim([0, 3])\nplt.xlabel('Frequency (Hz)')\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:44:38.835529Z","iopub.execute_input":"2022-10-26T21:44:38.835944Z","iopub.status.idle":"2022-10-26T21:44:39.20123Z","shell.execute_reply.started":"2022-10-26T21:44:38.835904Z","shell.execute_reply":"2022-10-26T21:44:39.199836Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n> To recover the hidden time information, Dennis Gabor introduced in the year 1946 the short-time Fourier transform (STFT). Instead of considering the entire signal, the main idea of the STFT is to consider only a small section of the signal. To this end, one fixes a so-called window function, which is a function that is nonzero for only a short period of time (defining the considered section). The original signal is then multiplied with the window function to yield a windowed signal. To obtain frequency information at different time instances, one shifts the window function across time and computes a Fourier transform for each of the resulting windowed signals. This idea is illustrated by the next example.","metadata":{}},{"cell_type":"code","source":"def windowed_ft(t, x, Fs, w_pos_sec, w_len):\n \n N = len(x)\n w_pos = int(Fs * w_pos_sec)\n w_padded = np.zeros(N)\n w_padded[w_pos:w_pos + w_len] = 1\n x = x * w_padded \n plt.figure(figsize=(8, 2))\n\n plt.subplot(1, 2, 1)\n plt.plot(t, x, c='k')\n plt.plot(t, w_padded, c='r')\n plt.xlim([min(t), max(t)])\n plt.ylim([-1.1, 1.1])\n plt.xlabel('Time (seconds)')\n\n plt.subplot(1, 2, 2)\n X = np.abs(np.fft.fft(x)) / Fs\n freq = np.fft.fftfreq(N, d=1/Fs)\n X = X[:N//2]\n freq = freq[:N//2]\n plt.plot(freq, X, c='k')\n plt.xlim([0, 7])\n plt.ylim([0, 3])\n plt.xlabel('Frequency (Hz)')\n plt.tight_layout()\n plt.show()\n \nw_len = 4 * Fs\nwindowed_ft(t, x, Fs, w_pos_sec=1, w_len=w_len)\nwindowed_ft(t, x, Fs, w_pos_sec=3, w_len=w_len)\nwindowed_ft(t, x, Fs, w_pos_sec=5, w_len=w_len)\n\nprint('Interactive interface for experimenting with different window shifts:')\ninteract(windowed_ft,\n w_pos_sec=FloatSlider(min=0, max=duration-(w_len/Fs), step=0.1, \n continuous_update=False, value=1.7, description='Position'),\n t=fixed(t), x=fixed(x), Fs=fixed(Fs), w_len=fixed(w_len));","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:44:39.202918Z","iopub.execute_input":"2022-10-26T21:44:39.203331Z","iopub.status.idle":"2022-10-26T21:44:40.48659Z","shell.execute_reply.started":"2022-10-26T21:44:39.203291Z","shell.execute_reply":"2022-10-26T21:44:40.485807Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> It is important to note that the STFT reflects not only the properties of the original signal but also those of the window function. First of all, the STFT depends on the length of the window, which determines the size of the section. Then, the STFT is influenced by the shape of the window. For example, the sharp edges of the rectangular window typically introduce \"ripple\" artifacts. We discuss such issues in more detail later.","metadata":{}},{"cell_type":"markdown","source":"
\n

Formal Definition of the Discrete STFT

We now consider the discrete case of the STFT and specify the most important mathematical formulas as needed in practical applications. Let x:[0:L&#x2212;1]:={0,1,&#x2026;,L&#x2212;1}&#x2192;R\" role=\"presentation\">x:[0:L−1]:={0,1,…,L−1}→Rx:[0:L−1]:={0,1,…,L−1}→R be a real-valued discrete-time (DT) signal of length L\" role=\"presentation\">LL obtained by equidistant sampling with respect to a fixed sampling rate Fs\" role=\"presentation\">FsFs given in Hertz. Furthermore, let w:[0:N&#x2212;1]&#x2192;R\" role=\"presentation\">w:[0:N−1]→Rw:[0:N−1]→R be a sampled window function of length N&#x2208;N\" role=\"presentation\">N∈NN∈N. For example, in the case of a rectangular window one has w(n)=1\" role=\"presentation\">w(n)=1w(n)=1 for n&#x2208;[0:N&#x2212;1]\" role=\"presentation\">n∈[0:N−1]n∈[0:N−1]. The length parameter N\" role=\"presentation\">NN determines the duration of the considered sections, which amounts to N/Fs\" role=\"presentation\">N/FsN/Fs seconds. One also introduces an additional parameter H&#x2208;N\" role=\"presentation\">H∈NH∈N, which is referred to as the hop size. The hop size parameter is specified in samples and determines the step size in which the window is to be shifted across the signal. With regard to these parameters, the discrete STFT X\" role=\"presentation\">XX of the signal x\" role=\"presentation\">xx is given by

\n
X(m,k):=&#x2211;n=0N&#x2212;1x(n+mH)w(n)exp(&#x2212;2&#x03C0;ikn/N)\" role=\"presentation\">X(m,k):=∑n=0N−1x(n+mH)w(n)exp(−2πikn/N)X(m,k):=∑n=0N−1x(n+mH)w(n)exp(−2πikn/N)

\n

with m&#x2208;[0:M]\" role=\"presentation\">m∈[0:M]m∈[0:M] and k&#x2208;[0:K]\" role=\"presentation\">k∈[0:K]k∈[0:K]. The number M:=&#x230A;L&#x2212;NH&#x230B;\" role=\"presentation\">M:=⌊L−NH⌋M:=⌊L−NH⌋ is the maximal frame index such that the window's time range is fully contained in the signal's time range. (We will see later some variants using padding strategies.) Furthermore, K=N/2\" role=\"presentation\">K=N/2K=N/2 (assuming that N\" role=\"presentation\">NN is even) is the frequency index corresponding to the Nyquist frequency. The complex number X(m,k)\" role=\"presentation\">X(m,k)X(m,k) denotes the kth\" role=\"presentation\">kthkth Fourier coefficient for the mth\" role=\"presentation\">mthmth time frame. Note that for each fixed time frame m\" role=\"presentation\">mm, one obtains a spectral vector of size K+1\" role=\"presentation\">K+1K+1 given by the coefficients X(m,k)\" role=\"presentation\">X(m,k)X(m,k) for k&#x2208;[0:K]\" role=\"presentation\">k∈[0:K]k∈[0:K]. The computation of each such spectral vector amounts to a DFT of size N\" role=\"presentation\">NN, which can be done efficiently using the FFT (see the FMP notebook on the FFT).

\n\n
","metadata":{}},{"cell_type":"markdown","source":"Spectrogram\nThe spectrogram is a two-dimensional representation of the squared magnitude of the STFT","metadata":{}},{"cell_type":"code","source":"def stft_basic(x, w, H=8, only_positive_frequencies=False):\n \"\"\"Compute a basic version of the discrete short-time Fourier transform (STFT)\n\n Notebook: C2/C2_STFT-Basic.ipynb\n\n Args:\n x (np.ndarray): Signal to be transformed\n w (np.ndarray): Window function\n H (int): Hopsize (Default value = 8)\n only_positive_frequencies (bool): Return only positive frequency part of spectrum (non-invertible)\n (Default value = False)\n\n Returns:\n X (np.ndarray): The discrete short-time Fourier transform\n \"\"\"\n N = len(w)\n L = len(x)\n M = np.floor((L - N) / H).astype(int) + 1\n X = np.zeros((N, M), dtype='complex')\n for m in range(M):\n x_win = x[m * H:m * H + N] * w\n X_win = np.fft.fft(x_win)\n X[:, m] = X_win\n\n if only_positive_frequencies:\n K = 1 + N // 2\n X = X[0:K, :]\n return X\n\nH = 8\nN = 128\nw = np.ones(N)\nX = stft_basic(x, w, H, only_positive_frequencies=True)\nY = np.abs(X) ** 2\n\nplt.figure(figsize=(8, 2))\nplt.subplot(1, 2, 1)\nplt.plot(np.arange(len(t)), x, c='k')\nplt.xlim([0, len(t)])\nplt.xlabel('Index (samples)')\nplt.subplot(1, 2, 2)\nplt.imshow(Y, origin='lower', aspect='auto', cmap='gray_r')\nplt.xlabel('Index (frames)')\nplt.ylabel('Index (frequency)')\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:44:40.48772Z","iopub.execute_input":"2022-10-26T21:44:40.488251Z","iopub.status.idle":"2022-10-26T21:44:40.860639Z","shell.execute_reply.started":"2022-10-26T21:44:40.488216Z","shell.execute_reply":"2022-10-26T21:44:40.859375Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n
\n
\n

Interpretation of Time and Frequency Indices

As for the temporal dimension, each Fourier coefficient X(m,k)\" role=\"presentation\">X(m,k)X(m,k) is associated with the physical time position

\n
Tcoef(m):=m&#x22C5;HFs\" role=\"presentation\">Tcoef(m):=m⋅HFsTcoef(m):=m⋅HFs

given in seconds. For example, for the smallest possible hop size H=1\" role=\"presentation\">H=1H=1, one obtains Tcoef(m)=m/Fs=m&#x22C5;T&#xA0;sec\" role=\"presentation\">Tcoef(m)=m/Fs=m⋅T secTcoef(m)=m/Fs=m⋅T sec. In this case, one obtains a spectral vector for each sample of the DT-signal x\" role=\"presentation\">xx, which results in a huge increase in data volume. Furthermore, considering sections that are only shifted by one sample generally yields very similar spectral vectors. To reduce this type of redundancy, one typically relates the hop size to the length N\" role=\"presentation\">NN of the window. For example, one often chooses H=N/2\" role=\"presentation\">H=N/2H=N/2, which constitutes a good trade-off between a reasonable temporal resolution and the data volume comprising all generated spectral coefficients. As for the frequency dimension, the index k\" role=\"presentation\">kk of X(m,k)\" role=\"presentation\">X(m,k)X(m,k) corresponds to the physical frequency

\n
Fcoef(k):=k&#x22C5;FsN\" role=\"presentation\">Fcoef(k):=k⋅FsNFcoef(k):=k⋅FsN

given in Hertz.

\n\n
\n
\n
\ngiven in Hertz.","metadata":{}},{"cell_type":"markdown","source":"\n\n","metadata":{}},{"cell_type":"code","source":"T_coef = np.arange(X.shape[1]) * H / Fs\nF_coef = np.arange(X.shape[0]) * Fs / N\n\nplt.figure(figsize=(8, 2))\n\nplt.subplot(1, 2, 1)\nplt.plot(t, x, c='k')\nplt.xlim([min(t), max(t)])\nplt.xlabel('Time (seconds)')\n\nplt.subplot(1, 2, 2)\nleft = min(T_coef)\nright = max(T_coef) + N / Fs\nlower = min(F_coef)\nupper = max(F_coef)\nplt.imshow(Y, origin='lower', aspect='auto', cmap='gray_r', \n extent=[left, right, lower, upper])\nplt.xlabel('Time (seconds)')\nplt.ylabel('Frequency (Hz)')\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:44:40.862553Z","iopub.execute_input":"2022-10-26T21:44:40.862901Z","iopub.status.idle":"2022-10-26T21:44:41.213291Z","shell.execute_reply.started":"2022-10-26T21:44:40.862869Z","shell.execute_reply":"2022-10-26T21:44:41.21223Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"https://www.audiolabs-erlangen.de/resources/MIR/FMP/C2/C2_STFT-Basic.html\nAcknowledgment: This Chapter was created by Meinard Müller and Frank Zalkow.\nhttps://www.audiolabs-erlangen.de/fau/professor/mueller\nhttps://www.audiolabs-erlangen.de/fau/assistant/zalkow","metadata":{}},{"cell_type":"markdown","source":"#### Audio File test","metadata":{}},{"cell_type":"code","source":"import matplotlib\n#fn_wav = os.path.join('..', 'data', 'C2', 'FMP_C2_F10.wav')\n\n#fn_wav =\"../input/piano-triads-wavset/piano_triads/A_min_4_0.wav\"\n#fn_wav =\"../input/emergency-vehicle-siren-sounds/sounds/ambulance/sound_106.wav\"\nfn_wav =\"../input/sample-wav-audio-files/Beethoven_Diabelli_Variation_No._13.wav\"\n\n\nx, Fs = librosa.load(fn_wav)\n\nH = 1024\nN = 2048\nw = np.hanning(N)\nX = stft_basic(x, w, H)\nY = np.abs(X) ** 2\neps = np.finfo(float).eps\nY_db = 10 * np.log10(Y + eps)\n\nT_coef = np.arange(X.shape[1]) * H / Fs\nF_coef = np.arange(X.shape[0]) * Fs / N\n\n\nfig = plt.figure(figsize=(8, 5))\n\ngs = matplotlib.gridspec.GridSpec(3, 2, height_ratios=[1, 2, 2], width_ratios=[100, 2])\nax1, ax2, ax3, ax4, ax5, ax6 = [plt.subplot(gs[i]) for i in range(6)]\n\nt = np.arange(len(x)) / Fs\nax1.plot(t, x, c='gray')\nax1.set_xlim([min(t), max(t)])\n\nax2.set_visible(False)\n\nleft = min(T_coef)\nright = max(T_coef) + N / Fs\nlower = min(F_coef)\nupper = max(F_coef)\n\nim1 = ax3.imshow(Y, origin='lower', aspect='auto', cmap='gray_r', \n extent=[left, right, lower, upper])\nim1.set_clim([0, 1000])\nax3.set_ylim([0, 5000])\nax3.set_ylabel('Frequency (Hz)')\ncbar = fig.colorbar(im1, cax=ax4)\nax4.set_ylabel('Magnitude (linear)', rotation=90)\n\nim2 = ax5.imshow(Y_db, origin='lower', aspect='auto', cmap='gray_r', \n extent=[left, right, lower, upper])\nim2.set_clim([-30, 20])\nax5.set_ylim([0, 5000])\nax5.set_xlabel('Time (seconds)')\nax5.set_ylabel('Frequency (Hz)')\ncbar = fig.colorbar(im2, cax=ax6)\nax6.set_ylabel('Magnitude (dB)', rotation=90)\n\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:44:41.214992Z","iopub.execute_input":"2022-10-26T21:44:41.21568Z","iopub.status.idle":"2022-10-26T21:44:44.414802Z","shell.execute_reply.started":"2022-10-26T21:44:41.215642Z","shell.execute_reply":"2022-10-26T21:44:44.413344Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Custom python STFT library __libfmp__\nThis implementation also includes padding options and realizes the centric view as used in the FMP notebooks. For details and similar options, we refer to the FMP notebook on conventions and implementations and the documentation of librosa. In the following code cell, we call the libfmp function.","metadata":{}},{"cell_type":"code","source":"pip install libfmp","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:44:44.416973Z","iopub.execute_input":"2022-10-26T21:44:44.418004Z","iopub.status.idle":"2022-10-26T21:45:26.803787Z","shell.execute_reply.started":"2022-10-26T21:44:44.417941Z","shell.execute_reply":"2022-10-26T21:45:26.802259Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import sys\nsys.path.append('..')\nimport libfmp.c2\n\n#fn_wav = os.path.join('..', 'data', 'C2', 'FMP_C2_F10.wav')\n#fn_wav = os.path.join('..', 'data', 'C2', 'FMP_C2_F10.wav')\n\n#fn_wav =\"../input/piano-triads-wavset/piano_triads/A_min_4_0.wav\"\nfn_wav =\"../input/sample-wav-audio-files/Mozart_from_Piano_Sonata_K310_first_movement.wav\"\n\n\nx, Fs = librosa.load(fn_wav)\n\nH = 1024\nN = 2048\nw = np.hanning(N)\nX = libfmp.c2.stft(x, w, H)\nY = np.abs(X) ** 2\neps = np.finfo(float).eps\nY_db = 10 * np.log10(Y + eps)\n\nT_coef = np.arange(X.shape[1]) * H / Fs\nF_coef = np.arange(X.shape[0]) * Fs / N\n\n\nfig = plt.figure(figsize=(8, 3))\n\ngs = matplotlib.gridspec.GridSpec(2, 2, height_ratios=[1, 2], width_ratios=[100, 2])\n\nax1, ax2, ax3, ax4 = [plt.subplot(gs[i]) for i in range(4)]\n\nt = np.arange(len(x)) / Fs\nax1.plot(t, x, c='gray')\nax1.set_xlim([min(t), max(t)])\n\nax2.set_visible(False)\n\nleft = min(T_coef)\nright = max(T_coef) + N / Fs\nlower = min(F_coef)\nupper = max(F_coef)\n\nim = ax3.imshow(Y_db, origin='lower', aspect='auto', cmap='gray_r', \n extent=[left, right, lower, upper])\nim.set_clim([-30, 20])\nax3.set_ylim([0, 5000])\nax3.set_xlabel('Time (seconds)')\nax3.set_ylabel('Frequency (Hz)')\ncbar = fig.colorbar(im, cax=ax4)\nax4.set_ylabel('Magnitude (dB)', rotation=90)\n\nplt.tight_layout()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:45:26.805518Z","iopub.execute_input":"2022-10-26T21:45:26.805937Z","iopub.status.idle":"2022-10-26T21:45:33.787953Z","shell.execute_reply.started":"2022-10-26T21:45:26.805892Z","shell.execute_reply":"2022-10-26T21:45:33.78659Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## [G2Net] Understand the Data\nby AYUSH THAKUR \n> https://www.kaggle.com/code/ayuraj/g2net-understand-the-data","metadata":{}},{"cell_type":"code","source":"!pip -qq install git+https://github.com/PyFstat/PyFstat@python37","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:45:33.789826Z","iopub.execute_input":"2022-10-26T21:45:33.791032Z","iopub.status.idle":"2022-10-26T21:45:51.792543Z","shell.execute_reply.started":"2022-10-26T21:45:33.790961Z","shell.execute_reply":"2022-10-26T21:45:51.790401Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# To create graitational waves\nimport pyfstat\nfrom pyfstat.utils import get_sft_as_arrays","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:45:51.795488Z","iopub.execute_input":"2022-10-26T21:45:51.796129Z","iopub.status.idle":"2022-10-26T21:45:51.80518Z","shell.execute_reply.started":"2022-10-26T21:45:51.796042Z","shell.execute_reply":"2022-10-26T21:45:51.803204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ROOT_DIR = '../input/g2net-detecting-continuous-gravitational-waves'\nos.path.isdir(ROOT_DIR)","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:45:51.806471Z","iopub.execute_input":"2022-10-26T21:45:51.806924Z","iopub.status.idle":"2022-10-26T21:45:51.92848Z","shell.execute_reply.started":"2022-10-26T21:45:51.806887Z","shell.execute_reply":"2022-10-26T21:45:51.926982Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.read_csv(f\"{ROOT_DIR}/train_labels.csv\")\ndf.head()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:45:51.93085Z","iopub.execute_input":"2022-10-26T21:45:51.932783Z","iopub.status.idle":"2022-10-26T21:45:51.97623Z","shell.execute_reply.started":"2022-10-26T21:45:51.932722Z","shell.execute_reply":"2022-10-26T21:45:51.974548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.target.value_counts()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:45:51.977992Z","iopub.execute_input":"2022-10-26T21:45:51.978747Z","iopub.status.idle":"2022-10-26T21:45:52.01072Z","shell.execute_reply.started":"2022-10-26T21:45:51.978697Z","shell.execute_reply":"2022-10-26T21:45:52.009447Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Read HDF5 File","metadata":{}},{"cell_type":"code","source":"from glob import glob\nfrom pathlib import Path\nimport h5py","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:45:52.012264Z","iopub.execute_input":"2022-10-26T21:45:52.012666Z","iopub.status.idle":"2022-10-26T21:45:52.06607Z","shell.execute_reply.started":"2022-10-26T21:45:52.01263Z","shell.execute_reply":"2022-10-26T21:45:52.064054Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_files = glob(f\"{ROOT_DIR}/train/*.hdf5\")\ntest_files = glob(f\"{ROOT_DIR}/test/*.hdf5\")\ntrain_files[:5]","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:45:52.068366Z","iopub.execute_input":"2022-10-26T21:45:52.069244Z","iopub.status.idle":"2022-10-26T21:45:52.151574Z","shell.execute_reply.started":"2022-10-26T21:45:52.0692Z","shell.execute_reply":"2022-10-26T21:45:52.150683Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DATA_IDX = 0","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:45:52.152924Z","iopub.execute_input":"2022-10-26T21:45:52.153552Z","iopub.status.idle":"2022-10-26T21:45:52.161655Z","shell.execute_reply.started":"2022-10-26T21:45:52.153515Z","shell.execute_reply":"2022-10-26T21:45:52.16019Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"file = Path(train_files[DATA_IDX])\nwith h5py.File(file, \"r\") as f:\n # Print all root level object names (aka keys) \n # these can be group or dataset names \n print(\"Top level Keys: %s\" % f.keys())\n key = list(f.keys())[0]\n print(\"This is the name of the file \\n\")\n # get the object type for a_group_key: usually group or dataset\n print(f\"Type of the data accessed using {key}: {type(f[key])}\")\n print(\"The type can be either Group or Dataset. Group has more keys while we can get numpy array from Dataset directly.\\n\")\n # If a_group_key is a group name, \n # this gets the object names in the group and returns as a list\n f = f[key]\n data_keys = list(f)\n print(\"These are the keys for the underlying group: \", data_keys, \"\\n\")\n # get type of each data_keys\n print(type(f[data_keys[0]])) # Group\n print(type(f[data_keys[1]])) # Group\n print(type(f[data_keys[2]])) # Dataset\n print(\"The underlying data type for each key are Group, Group and Dataset. We thus need to go one level deep to access the data from Group. \\n\")\n # H1\n h1 = f[data_keys[0]]\n h1_keys = list(h1)\n print(\"These are the keys for the Group accessed uisng H1 key: \", h1_keys, \"\\n\")\n # L1\n l1 = f[data_keys[1]]\n l1_keys = list(l1)\n print(\"These are the keys for the Group accessed uisng L1 key: \", l1_keys, \"\\n\")\n # frequency_Hz\n freq_hz = f[data_keys[2]]\n freq_hz = list(freq_hz)\n print(\"Since the type accessed using frequency_hz key was Dataset, we can get the array directly. The length of the array is: \", len(freq_hz), \"\\n\")\n # H1 data\n h1_stft = h1[\"SFTs\"][()]\n h1_timestamp = h1[\"timestamps_GPS\"][()]\n print(\"The H1 SFTs data shape: \", h1_stft.shape)\n print(\"The H1 timestamp data shape: \", h1_timestamp.shape, \"\\n\")\n # H2 data\n l1_stft = l1[\"SFTs\"][()]\n l1_timestamp = l1[\"timestamps_GPS\"][()]\n print(\"The L1 SFTs data shape: \", l1_stft.shape)\n print(\"The L1 timestamp data shape: \", l1_timestamp.shape, \"\\n\")\n print(\"We have time on the x-axis while frequency on the y-axis. \\n\")\n \n print(\"Label: \", df.loc[DATA_IDX].target)","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:45:52.163619Z","iopub.execute_input":"2022-10-26T21:45:52.164145Z","iopub.status.idle":"2022-10-26T21:45:52.677344Z","shell.execute_reply.started":"2022-10-26T21:45:52.164073Z","shell.execute_reply":"2022-10-26T21:45:52.676034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Utility to read hdf5 file","metadata":{}},{"cell_type":"code","source":"# Utility to read hdf5 file\ndef read_data(file: Path):\n with h5py.File(file, \"r\") as f:\n filename = file.stem\n f = f[filename]\n h1 = f[\"H1\"]\n l1 = f[\"L1\"]\n freq_hz = list(f[\"frequency_Hz\"])\n \n h1_stft = h1[\"SFTs\"][()]\n h1_timestamp = h1[\"timestamps_GPS\"][()]\n # H2 data\n l1_stft = l1[\"SFTs\"][()]\n l1_timestamp = l1[\"timestamps_GPS\"][()]\n \n return {\n \"H1\": [h1_stft, h1_timestamp],\n \"L1\": [l1_stft, l1_timestamp],\n \"freq_hz\": freq_hz\n }","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:45:52.686849Z","iopub.execute_input":"2022-10-26T21:45:52.688038Z","iopub.status.idle":"2022-10-26T21:45:52.700931Z","shell.execute_reply.started":"2022-10-26T21:45:52.687985Z","shell.execute_reply":"2022-10-26T21:45:52.694976Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Shapes","metadata":{}},{"cell_type":"code","source":"filename=Path(train_files[1])\ndata = read_data(filename)\nh1_sft, h1_ts = data[\"H1\"]\nl1_sft, l1_ts = data[\"L1\"]\nfreq_hz = data[\"freq_hz\"]\nprint(h1_sft.shape, h1_ts.shape, len(freq_hz))","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:45:53.232362Z","iopub.execute_input":"2022-10-26T21:45:53.232848Z","iopub.status.idle":"2022-10-26T21:45:53.426869Z","shell.execute_reply.started":"2022-10-26T21:45:53.232795Z","shell.execute_reply":"2022-10-26T21:45:53.4252Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_shapes=pd.read_csv(\"../input/g2net-datashape/train_features.csv\")\ndf_shapes.tail()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:45:53.436241Z","iopub.execute_input":"2022-10-26T21:45:53.436911Z","iopub.status.idle":"2022-10-26T21:45:53.473628Z","shell.execute_reply.started":"2022-10-26T21:45:53.436874Z","shell.execute_reply":"2022-10-26T21:45:53.472677Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"df_shapes.describe()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:45:53.475445Z","iopub.execute_input":"2022-10-26T21:45:53.476303Z","iopub.status.idle":"2022-10-26T21:45:53.514383Z","shell.execute_reply.started":"2022-10-26T21:45:53.476259Z","shell.execute_reply":"2022-10-26T21:45:53.513163Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":" #### Outliers in trainset\nThere are some data points that differs significantly from other observations. An outlier may be due to variability in the measurement or it may indicate experimental error; the latter are sometimes excluded from the data set. An outlier can cause serious problems in statistical analyses.","metadata":{}},{"cell_type":"markdown","source":" ###### LIGO Hanford Observatory observations","metadata":{}},{"cell_type":"code","source":"sns.displot(df_shapes, x=\"h1_ts\", bins=20)\nplt.figure(figsize=(4,3))","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:45:53.516139Z","iopub.execute_input":"2022-10-26T21:45:53.516501Z","iopub.status.idle":"2022-10-26T21:45:53.954056Z","shell.execute_reply.started":"2022-10-26T21:45:53.516467Z","shell.execute_reply":"2022-10-26T21:45:53.952685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"###### LIGO Livingston","metadata":{}},{"cell_type":"code","source":"df_shapeL=pd.read_csv(\"../input/g2net-datashape/Ltrain_features.csv\")\n","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:20.257734Z","iopub.execute_input":"2022-10-26T21:46:20.258265Z","iopub.status.idle":"2022-10-26T21:46:20.26858Z","shell.execute_reply.started":"2022-10-26T21:46:20.258223Z","shell.execute_reply":"2022-10-26T21:46:20.267147Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.displot(df_shapeL, x=\"h1_ts\", bins=20)\nplt.figure(figsize=(4,3))","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:21.458602Z","iopub.execute_input":"2022-10-26T21:46:21.459127Z","iopub.status.idle":"2022-10-26T21:46:21.825162Z","shell.execute_reply.started":"2022-10-26T21:46:21.459068Z","shell.execute_reply":"2022-10-26T21:46:21.823991Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Outliers in testset","metadata":{}},{"cell_type":"code","source":"df_shapes_test=pd.read_csv(\"../input/g2net-datashape/test_features.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:34.087963Z","iopub.execute_input":"2022-10-26T21:46:34.088492Z","iopub.status.idle":"2022-10-26T21:46:34.131792Z","shell.execute_reply.started":"2022-10-26T21:46:34.088444Z","shell.execute_reply":"2022-10-26T21:46:34.130563Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.displot(df_shapes_test, x=\"h1_ts\", bins=20)\nplt.figure(figsize=(4,3))","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:35.390313Z","iopub.execute_input":"2022-10-26T21:46:35.391455Z","iopub.status.idle":"2022-10-26T21:46:35.776084Z","shell.execute_reply.started":"2022-10-26T21:46:35.391408Z","shell.execute_reply":"2022-10-26T21:46:35.775021Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.distplot(df_shapes_test[['h1_ts']], hist=False, rug=True)\nsns.distplot(df_shapes[['h1_ts']], hist=False, rug=True)\nsns.distplot(df_shapesL[['h1_ts']], hist=False, rug=True)\n\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:35.777961Z","iopub.execute_input":"2022-10-26T21:46:35.77856Z","iopub.status.idle":"2022-10-26T21:46:36.18502Z","shell.execute_reply.started":"2022-10-26T21:46:35.77852Z","shell.execute_reply":"2022-10-26T21:46:36.184142Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.distplot(df_shapes_test[['h1_ts']], hist=False, rug=True).set_xlim(4200,4900)\nsns.distplot(df_shapes[['h1_ts']], hist=False, rug=True).set_xlim(4200,4900)\nsns.distplot(df_shapesL[['h1_ts']], hist=False, rug=True).set_xlim(4200,4900)\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:36.186644Z","iopub.execute_input":"2022-10-26T21:46:36.18721Z","iopub.status.idle":"2022-10-26T21:46:36.747423Z","shell.execute_reply.started":"2022-10-26T21:46:36.187174Z","shell.execute_reply":"2022-10-26T21:46:36.746144Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_shapes_test.to_csv(\"test_features.csv\")\n#df_shapes.to_csv(\"train_features.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:36.749632Z","iopub.execute_input":"2022-10-26T21:46:36.750288Z","iopub.status.idle":"2022-10-26T21:46:36.790578Z","shell.execute_reply.started":"2022-10-26T21:46:36.750242Z","shell.execute_reply":"2022-10-26T21:46:36.788903Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del df_shapes_test\ndel df_shapes","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:37.3732Z","iopub.execute_input":"2022-10-26T21:46:37.374276Z","iopub.status.idle":"2022-10-26T21:46:37.380537Z","shell.execute_reply.started":"2022-10-26T21:46:37.374221Z","shell.execute_reply":"2022-10-26T21:46:37.379184Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### H1: LIGO Hanford","metadata":{}},{"cell_type":"code","source":"from datetime import datetime\nprint(datetime.utcfromtimestamp(h1_ts[0]).strftime('%Y-%m-%d %H:%M:%S'))","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:38.040893Z","iopub.execute_input":"2022-10-26T21:46:38.042512Z","iopub.status.idle":"2022-10-26T21:46:38.050761Z","shell.execute_reply.started":"2022-10-26T21:46:38.042454Z","shell.execute_reply":"2022-10-26T21:46:38.049105Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### STFT LIGO Hanford\nShort time fourier transform computes the fourier transform of the signal at shorter duration of signals given by a window size. The DFT (discrete fourier transform) or FFT (fast fourier transform) are computed by taking the entire duration of the signal into consideration.\n\nNote that the plots generated below are not considering the gaps in the timestamp and considers the data as a single continuous SFT","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots(2, 2, figsize=(16, 10))\nfig.suptitle(f\"Filename Id: {file.stem}\")\n\nfor d_ind, detector in enumerate([\"H1\", \"L1\"]):\n ax[d_ind][0].set(xlabel=\"Timestamps [GPS]\",\n ylabel=\"Frequency [Hz]\",\n title=f\"{detector} - Real part\")\n ax[d_ind][1].set(xlabel=\"Timestamps [GPS]\",\n ylabel=\"Frequency [Hz]\",\n title=f\"{detector} - Imaginary part\")\n\n c0 = ax[d_ind][0].pcolormesh(data[detector][1], data[\"freq_hz\"], \n data[detector][0].real)\n c1 = ax[d_ind][1].pcolormesh(data[detector][1], data[\"freq_hz\"], \n data[detector][0].imag)\n\n fig.colorbar(c0, ax=ax[d_ind][0])\n fig.colorbar(c1, ax=ax[d_ind][1])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:38.949694Z","iopub.execute_input":"2022-10-26T21:46:38.950176Z","iopub.status.idle":"2022-10-26T21:46:43.588651Z","shell.execute_reply.started":"2022-10-26T21:46:38.950136Z","shell.execute_reply":"2022-10-26T21:46:43.587475Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Understanding By Generating","metadata":{}},{"cell_type":"markdown","source":"The data can better be understood if we generate the data and see how different parameters contribute to the SFTs. This section is based on the awesome tutorials provided by the host of this competition. You can find the tutorials here. The following is a distilled information in an attemp to better understand the data.","metadata":{}},{"cell_type":"markdown","source":"Understanding By Generating\nThe data can better be understood if we generate the data and see how different parameters contribute to the SFTs. This section is based on the awesome tutorials provided by the host of this competition. You can find the tutorials here. The following is a distilled information in an attemp to better understand the data.\n\nTLDR;\nThe SFTs generated with h0 parameter equal to 0 have target/label as 0. The SFTs generated with h0 parameter equal to 1, thus having SNR (signal to noise ratio) > 0 have it's target/label as 1. The competition hosts were kind enough to clarify this here.\n\nThe main piece of code to generate the data is writer = pyfstat.Writer(**writer_kwargs, **signal_parameters). If you don't pass **signal_parameter you will get data corresponding to 0 and so on. You can check the Writer class here.","metadata":{}},{"cell_type":"markdown","source":"#### CW wave properties\nThe CW wave has the following properties:\n\n__Quasi-monochromatic__: \"Wave in which most of the energy is confined to a single wavelength or a very narrow waveband\". However, gravitational continuous waves are not completely monochromatic. (source)\n**Long-standing**: The standing waves are stationary (not a travelling wave). The disturbance does not travel in any direction. Standing waves have points of zero amplitude called nodes and points of maximum amplitude called the antinodes.\n**Produced by non-axisymmetric rapidly-spinning neutron stars**\n**Produce narrow-banded signals** - that's why the y axis (Freq in Hz) in the spectrogram lies in **a narrow band of 1-2 Hz**.\nThe waves are **doppler-modulated** due to the motion of the detector in the Solar system.\nIt's a common practice to use short time fourier transform to analyze these waves. These are Fourier transforms of short data segments (**typically around 30 minutes or 1800 seconds**).","metadata":{}},{"cell_type":"markdown","source":"#### Generate Gausian Noise","metadata":{}},{"cell_type":"code","source":"# Setup Writer\nwriter_kwargs = {\n \"label\": \"single_detector_gaussian_noise\",\n \"outdir\": \"PyFstat_example_data\",\n \"tstart\": 1238166018, # Starting time of the observation [GPS time]\n \"duration\": 5 * 86400, # Duration [seconds]\n \"detectors\": \"H1\", # Detector to simulate, in this case LIGO Hanford\n \"F0\": 100.0, # Central frequency of the band to be generated [Hz]\n \"Band\": 1.0, # Frequency band-width around F0 [Hz]\n \"sqrtSX\": 1e-23, # Single-sided Amplitude Spectral Density of the noise\n \"Tsft\": 1800, # Fourier transform time duration\n \"SFTWindowType\": \"tukey\", # Window function to compute short Fourier transforms\n \"SFTWindowBeta\": 0.01, # Parameter associated to the window function\n}\nwriter = pyfstat.Writer(**writer_kwargs)\n\n# Create SFTs\nwriter.make_data()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:43.590638Z","iopub.execute_input":"2022-10-26T21:46:43.59108Z","iopub.status.idle":"2022-10-26T21:46:45.221597Z","shell.execute_reply.started":"2022-10-26T21:46:43.591047Z","shell.execute_reply":"2022-10-26T21:46:45.219857Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Note:** \"*As a side note, the path where the SFTs are created is stored in the sftfilepath attribute of Writer **(mind that it gets overwritten every time make_data is called!)**.\" (source)*","metadata":{}},{"cell_type":"code","source":"frequency, timestamps, fourier_data = get_sft_as_arrays(writer.sftfilepath)","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:45.223696Z","iopub.execute_input":"2022-10-26T21:46:45.224179Z","iopub.status.idle":"2022-10-26T21:46:45.249692Z","shell.execute_reply.started":"2022-10-26T21:46:45.224125Z","shell.execute_reply":"2022-10-26T21:46:45.248366Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"frequency, timestamps, fourier_data = get_sft_as_arrays(writer.sftfilepath)","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:45.252981Z","iopub.execute_input":"2022-10-26T21:46:45.254394Z","iopub.status.idle":"2022-10-26T21:46:45.274983Z","shell.execute_reply.started":"2022-10-26T21:46:45.254348Z","shell.execute_reply":"2022-10-26T21:46:45.273444Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The shape of H1 SFT is: (1800, 240) while that of it's time stamp is: (240,)\nThe shape of frequency Hz is: (1800,)\nWhy the shape of time stamp is (240,)?\n\nThere are three time related parameters:\n> * \"tstart\": 1238166018, # Starting time of the observation [GPS time]\n> * \"duration\": 5 * 86400, # Duration [seconds]\n> * \"Tsft\": 1800, # Fourier transform time duration","metadata":{}},{"cell_type":"code","source":"print(f\"Total duration is : {5*86400}; time segments for SFT computation: {1800}; total timestamps: {(5*86400)/1800}\")","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:45.276934Z","iopub.execute_input":"2022-10-26T21:46:45.277374Z","iopub.status.idle":"2022-10-26T21:46:45.286241Z","shell.execute_reply.started":"2022-10-26T21:46:45.277334Z","shell.execute_reply":"2022-10-26T21:46:45.284813Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\"\"\"\nUtils\n=====\n\nUtility functions to simplify tutorials.\n\"\"\"\nimport matplotlib.pyplot as plt\nimport numpy as np\nfrom matplotlib import colors\nfrom scipy import stats\n\nplt.rcParams[\"font.family\"] = \"serif\"\nplt.rcParams[\"font.size\"] = 20\n\n\ndef plot_real_imag_spectrograms(timestamps, frequency, fourier_data):\n \"uses simple plt.pcolormesh\"\n fig, axs = plt.subplots(1, 2, figsize=(16, 10))\n\n for ax in axs:\n ax.set(xlabel=\"SFT index\", ylabel=\"Frequency [Hz]\")\n\n time_in_days = (timestamps - timestamps[0]) / 1800\n\n axs[0].set_title(\"SFT Real part\")\n c = axs[0].pcolormesh(\n time_in_days,\n frequency,\n fourier_data.real,\n norm=colors.CenteredNorm(),\n )\n fig.colorbar(c, ax=axs[0], orientation=\"horizontal\", label=\"Fourier Amplitude\")\n\n axs[1].set_title(\"SFT Imaginary part\")\n c = axs[1].pcolormesh(\n time_in_days,\n frequency,\n fourier_data.imag,\n norm=colors.CenteredNorm(),\n )\n\n fig.colorbar(c, ax=axs[1], orientation=\"horizontal\", label=\"Fourier Amplitude\")\n\n return fig, axs\n\n\ndef plot_real_imag_spectrograms_with_gaps(timestamps, frequency, fourier_data, Tsft):\n\n # Fill up gaps with Nans\n gap_length = timestamps[1:] - (timestamps[:-1] + Tsft)\n\n gap_data = [fourier_data[:, 0]]\n gap_timestamps = [timestamps[0]]\n\n for ind, gap in enumerate(gap_length):\n if gap > 0:\n gap_data.append(np.full_like(fourier_data[:, ind], np.nan + 1j * np.nan))\n gap_timestamps.append(timestamps[ind] + Tsft)\n\n gap_data.append(fourier_data[:, ind + 1])\n gap_timestamps.append(timestamps[ind + 1])\n\n return plot_real_imag_spectrograms(\n np.hstack(gap_timestamps), frequency, np.vstack(gap_data).T\n )\n\n\ndef plot_real_imag_histogram(fourier_data, theoretical_stdev=None):\n\n fig, ax = plt.subplots(figsize=(16, 10))\n ax.set(xlabel=\"SFT value\", ylabel=\"PDF\", yscale=\"log\")\n\n ax.hist(\n fourier_data.real.ravel(),\n density=True,\n bins=\"auto\",\n histtype=\"step\",\n lw=2,\n label=\"Real part\",\n )\n ax.hist(\n fourier_data.imag.ravel(),\n density=True,\n bins=\"auto\",\n histtype=\"step\",\n lw=2,\n label=\"Imaginary part\",\n )\n\n if theoretical_stdev is not None:\n x = np.linspace(-4 * theoretical_stdev, 4 * theoretical_stdev, 1000)\n y = stats.norm(scale=theoretical_stdev).pdf(x)\n ax.plot(x, y, color=\"black\", ls=\"--\", label=\"Gaussian distribution\")\n\n ax.legend()\n\n return fig, ax\n\n\ndef plot_amplitude_phase_spectrograms(timestamps, frequency, fourier_data):\n fig, axs = plt.subplots(1, 2, figsize=(16, 10))\n\n for ax in axs:\n ax.set(xlabel=\"SFT index\", ylabel=\"Frequency [Hz]\")\n\n time_in_days = (timestamps - timestamps[0]) / 1800\n\n axs[0].set_title(\"SFT absolute value\")\n c = axs[0].pcolorfast(\n time_in_days, frequency, np.absolute(fourier_data), norm=colors.Normalize()\n )\n fig.colorbar(c, ax=axs[0], orientation=\"horizontal\", label=\"Value\")\n\n axs[1].set_title(\"SFT phase\")\n c = axs[1].pcolorfast(\n time_in_days, frequency, np.angle(fourier_data), norm=colors.CenteredNorm()\n )\n\n fig.colorbar(c, ax=axs[1], orientation=\"horizontal\", label=\"Value\")\n\n return fig, axs","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:45.288766Z","iopub.execute_input":"2022-10-26T21:46:45.28937Z","iopub.status.idle":"2022-10-26T21:46:45.319684Z","shell.execute_reply.started":"2022-10-26T21:46:45.289315Z","shell.execute_reply":"2022-10-26T21:46:45.31823Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_real_imag_spectrograms(\n timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n);","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:45.321407Z","iopub.execute_input":"2022-10-26T21:46:45.322065Z","iopub.status.idle":"2022-10-26T21:46:46.931649Z","shell.execute_reply.started":"2022-10-26T21:46:45.322028Z","shell.execute_reply":"2022-10-26T21:46:46.930201Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**The SNR for such a signal is 0 since there are no signal. :P :P :)**\n\nLet's inspect the apmplitude and angle of the generate SFT.\nSince it's a gaussian noise, the power (amplitude) and angle in the SFT is random.","metadata":{}},{"cell_type":"code","source":"fig, ax = plot_amplitude_phase_spectrograms(\n timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n);","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:46.933256Z","iopub.execute_input":"2022-10-26T21:46:46.934113Z","iopub.status.idle":"2022-10-26T21:46:47.705929Z","shell.execute_reply.started":"2022-10-26T21:46:46.934049Z","shell.execute_reply":"2022-10-26T21:46:47.704947Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#https://kevinsprojects.wordpress.com/2014/12/13/short-time-fourier-transform-using-python-and-numpy/","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:47.707168Z","iopub.execute_input":"2022-10-26T21:46:47.70823Z","iopub.status.idle":"2022-10-26T21:46:47.712466Z","shell.execute_reply.started":"2022-10-26T21:46:47.708189Z","shell.execute_reply":"2022-10-26T21:46:47.711429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Generate CW Signal\n​\nThe CW waves are primarily dependent on: `amplitude` and `Doppler` parameters.\n​\n- the `amplitide` parameter depends on nominal CW amplitude $h0$, the cosine of the source inclination angle with respect to the line of sight $cos(ι)$ (cosi), the polarization angle $\\psi{}$ (psi) and the initial phase of the wave $\\phi{}$ (phi). Depending on the emission mechanism, $h0$ can be further described using further physical quantities such as the source's frequency, ellipticity, or distance to the detector.\n​\n- the `doppler` parameter describe the relationship between the source of emmission (and it's motion) and the detectors. The parameter depends on gravitational wave frequency $f0$, spindown $f1$and sky position $\\hat{n}$, which we will parametrize using the right ascension $\\alpha{}$ and declination $\\delta{}$ angles.","metadata":{}},{"cell_type":"markdown","source":"**In this case the signal is really clear and thus SNR should be high.**","metadata":{}},{"cell_type":"code","source":"# The same writer params are required to generate gaussian noise.\nwriter_kwargs = {\n \"label\": \"single_detector_gaussian_noise\",\n \"outdir\": \"PyFstat_example_data\",\n \"tstart\": 1238166018,\n \"duration\": 365 * 86400,\n \"detectors\": \"H1\",\n \"sqrtSX\": 1e-23,\n \"Tsft\": 1800,\n \"SFTWindowType\": \"tukey\",\n \"SFTWindowBeta\": 0.01,\n}\n\n# These parameters determine the CW signal that we have to detect.\nsignal_parameters = {\n \"F0\": 100.0,\n \"F1\": -1e-9,\n \"Alpha\": 0.0,\n \"Delta\": 0.0,\n \"h0\": 1e-22, # this parameter in particular determines if a signal exist in the generated SFT.\n \"cosi\": 1,\n \"psi\": 0.0,\n \"phi\": 0.0,\n \"tref\": writer_kwargs[\"tstart\"],\n}\n\nwriter = pyfstat.Writer(**writer_kwargs, **signal_parameters)\nwriter.make_data()\nfrequency, timestamps, fourier_data = get_sft_as_arrays(writer.sftfilepath)\nplot_real_imag_spectrograms(\n timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n);","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:46:47.715831Z","iopub.execute_input":"2022-10-26T21:46:47.716543Z","iopub.status.idle":"2022-10-26T21:47:01.480479Z","shell.execute_reply.started":"2022-10-26T21:46:47.716504Z","shell.execute_reply":"2022-10-26T21:47:01.479186Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The power and the angle component of the generate SFT is prominent.","metadata":{}},{"cell_type":"markdown","source":"Change the variables","metadata":{}},{"cell_type":"code","source":"# The same writer params are required to generate gaussian noise.\nwriter_kwargs = {\n \"label\": \"single_detector_gaussian_noise\",\n \"outdir\": \"PyFstat_example_data\",\n \"tstart\": 1238166018,\n \"duration\": 365 * 86400,\n \"detectors\": \"H1\",\n \"sqrtSX\": 1e-23,\n \"Tsft\": 1800,\n \"SFTWindowType\": \"tukey\",\n \"SFTWindowBeta\": 0.01,\n}\n\n# These parameters determine the CW signal that we have to detect.\nsignal_parameters = {\n \"F0\": 500.0,\n \"F1\": -1e-27,\n \"Alpha\": 3.,\n \"Delta\": 0.5,\n \"h0\": 1e-22, # this parameter in particular determines if a signal exist in the generated SFT. < \"sqrtSX\"\n \"cosi\": 0.1,\n \"psi\": 0.0,\n \"phi\": 2.,\n \"tref\": writer_kwargs[\"tstart\"],\n}\n\nwriter = pyfstat.Writer(**writer_kwargs, **signal_parameters)\nwriter.make_data()\nfrequency, timestamps, fourier_data = get_sft_as_arrays(writer.sftfilepath)\nplot_real_imag_spectrograms(\n timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n);","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:47:01.482679Z","iopub.execute_input":"2022-10-26T21:47:01.483663Z","iopub.status.idle":"2022-10-26T21:47:18.314813Z","shell.execute_reply.started":"2022-10-26T21:47:01.483606Z","shell.execute_reply":"2022-10-26T21:47:18.3138Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# SNR can be compute from a set of SFTs for a specific set\n# of parameters as follows:\nsnr = pyfstat.SignalToNoiseRatio.from_sfts(\n F0=writer.F0, sftfilepath=writer.sftfilepath\n)\nsquared_snr = snr.compute_snr2(\n Alpha=writer.Alpha, \n Delta=writer.Delta,\n psi=writer.psi,\n phi=writer.phi, \n h0=writer.h0,\n cosi=writer.cosi\n)\nsqrt_snr = np.sqrt(squared_snr)\nprint(f\"Signal - SNR: {sqrt_snr}\")","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:47:18.316498Z","iopub.execute_input":"2022-10-26T21:47:18.317437Z","iopub.status.idle":"2022-10-26T21:47:20.505267Z","shell.execute_reply.started":"2022-10-26T21:47:18.31739Z","shell.execute_reply":"2022-10-26T21:47:20.503877Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> The data of interest (real data) should have such high SNRs (there would be no point of this competition, if signal is easily detected).","metadata":{}},{"cell_type":"markdown","source":"### **SNR range of CW signal**\nLet's talk about the possible SNR range of CW signals. As per the tutorial shared by competition hosts, \"The presence of visible CW features in the spectrogram can be quantified using the so-called sensitivity depth, which is defined as the ratio between the noise and CW amplitude $D = \\sqrt{Sn}/h0$. Current searches for CW signals from unknown sources are able to detect a sensitivity depth ranging between 10 $Hz^{-0.5}$ and 50 $Hz^{-0.5}$. Visual features, however, tend to disappear around a depth of 20 $Hz^{-0.5}$.\"\n​\nLooking at the formula, $\\sqrt{Sn}$ is \"sqrtSX\" which in our case is `1e-23`. And we want to look at visual features where the sensitivity depth is 20 $Hz^{-0.5}$. Thus we will calculate $h0$.","metadata":{}},{"cell_type":"code","source":"# CW signal at a depth of 20\nsignal_parameters[\"h0\"] = writer_kwargs[\"sqrtSX\"] / 20.0\n\nwriter = pyfstat.Writer(**writer_kwargs, **signal_parameters)\n\n# Create SFTs\nwriter.make_data()\n\nfrequency, timestamps, fourier_data = get_sft_as_arrays(writer.sftfilepath)\n\nplot_real_imag_spectrograms(\n timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n)\n\nplot_amplitude_phase_spectrograms(\n timestamps[\"H1\"], frequency, fourier_data[\"H1\"]\n);","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:47:20.506975Z","iopub.execute_input":"2022-10-26T21:47:20.507534Z","iopub.status.idle":"2022-10-26T21:47:39.138714Z","shell.execute_reply.started":"2022-10-26T21:47:20.507482Z","shell.execute_reply":"2022-10-26T21:47:39.137156Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> We really can't see any visual feature in the spectrogram.","metadata":{}},{"cell_type":"code","source":"# SNR can be compute from a set of SFTs for a specific set\n# of parameters as follows:\nsnr = pyfstat.SignalToNoiseRatio.from_sfts(\n F0=writer.F0, sftfilepath=writer.sftfilepath\n)\nsquared_snr = snr.compute_snr2(\n Alpha=writer.Alpha,\n Delta=writer.Delta,\n psi=writer.psi,\n phi=writer.phi,\n h0=writer.h0,\n cosi=writer.cosi\n)\nsqrt_snr = np.sqrt(squared_snr)\nprint(f\"Signal - SNR: {sqrt_snr}\")","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:47:39.140709Z","iopub.execute_input":"2022-10-26T21:47:39.14173Z","iopub.status.idle":"2022-10-26T21:47:41.342607Z","shell.execute_reply.started":"2022-10-26T21:47:39.14168Z","shell.execute_reply":"2022-10-26T21:47:41.341208Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**But the SNR > 0 indicates that our target signal (CW waves) is present and we need better system to predict them.**","metadata":{}},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"## Imports 🗂","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%capture\n!pip install nexusformat","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:47:41.348595Z","iopub.execute_input":"2022-10-26T21:47:41.349206Z","iopub.status.idle":"2022-10-26T21:47:54.88574Z","shell.execute_reply.started":"2022-10-26T21:47:41.349149Z","shell.execute_reply":"2022-10-26T21:47:54.883842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\nfrom pathlib import Path\nimport os\nimport random\nimport numpy as np\nfrom pprint import pprint\nfrom tqdm.notebook import tqdm\n\nimport seaborn as sns\nsns.set_theme()\n\nfrom scipy import signal\nfrom scipy.fft import fftshift\nimport matplotlib.pyplot as plt\n\nimport h5py\nimport nexusformat.nexus as nx\n\nimport warnings\nwarnings.filterwarnings('ignore')\n\nfrom IPython.display import HTML, display\ndisplay(HTML(''))","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:47:54.88797Z","iopub.execute_input":"2022-10-26T21:47:54.888447Z","iopub.status.idle":"2022-10-26T21:47:54.941813Z","shell.execute_reply.started":"2022-10-26T21:47:54.888397Z","shell.execute_reply":"2022-10-26T21:47:54.940585Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"DATA_PATH = Path('../input/g2net-detecting-continuous-gravitational-waves')\nTRAIN_PATH = DATA_PATH/'train'\nTEST_PATH = DATA_PATH/'test'\ntrain_example_with_signal_path = TRAIN_PATH/'cc561e4fc.hdf5'\ntrain_example_without_signal_path = TRAIN_PATH/'fb6db0d08.hdf5'","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:47:54.943361Z","iopub.execute_input":"2022-10-26T21:47:54.943759Z","iopub.status.idle":"2022-10-26T21:47:54.950838Z","shell.execute_reply.started":"2022-10-26T21:47:54.943716Z","shell.execute_reply":"2022-10-26T21:47:54.949396Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# HDF5 🤔","metadata":{}},{"cell_type":"markdown","source":" The data is stored in HDF5 files. Hierarchical Data Format version 5 (HDF5), is an open source file format that supports large, complex, heterogeneous data. HDF5 uses a \"file directory\" like structure that allows you to organize data.\n\n HDF5 simplifies the file structure to include only two major types of object:\n- Datasets, which are typed multidimensional arrays\n- Groups, which are container structures that can hold datasets and other groups\n\n We can view the structure of our data with the nexusformat package.\n\n The `h5py` package is a thin, pythonic wrapper around HDF5 we can use it to quickly load our data.","metadata":{}},{"cell_type":"markdown","source":"
\"Space\"","metadata":{}},{"cell_type":"markdown","source":"## We can see the structure of the training and test data below:","metadata":{}},{"cell_type":"markdown","source":"- `ID` is the top group of the HDF5 file and links the datapoint to it's label in the `train_labels` csv (group)\n\n- `frequency_Hz` contains the range frequencies measured by the dectors (dataset)\n\n\n- `H1` contains the data for the LIGO Hanford decector (group) \n \n - `SFTs` is the Short-time Fourier Transforms amplitudes for each timestamp at each frequency (dataset)\n - `timestamps` contains the timestamps for the measurement (dataset)\n\n \n- `L1` contains the data for the LIGO Livingston decector (group) \n \n - `SFTs` is the Short-time Fourier Transforms amplitudes for each timestamp at each frequency (dataset)\n - `timestamps` contains the timestamps for the measurement (dataset) \n \n This structure can be visualised below 👇 ","metadata":{}},{"cell_type":"markdown","source":"![](https://i.imgur.com/M6xfOri.png)","metadata":{}},{"cell_type":"markdown","source":"## Let's look at an example","metadata":{}},{"cell_type":"code","source":"with h5py.File(train_example_with_signal_path, \"r\") as f:\n \n # get first object name/key; this is the data point ID\n ID_key = list(f.keys())[0]\n print(f\"ID: {ID_key} \\n\")\n \n # Retrieve the Livingston decector data\n print(f\"- {list(f[ID_key].keys())[1]}\")\n L1_SFTs = f[ID_key]['L1']['SFTs']\n print(f\"-- SFTs amplitudes: {L1_SFTs.shape}\")\n L1_ts = f[ID_key]['L1']['timestamps_GPS']\n print(f\"-- timestamps: {L1_ts.shape} \\n\")\n \n # Retrieve the Hanford decector data\n print(f\"- {list(f[ID_key].keys())[0]}\")\n H1_SFTs = f[ID_key]['H1']['SFTs']\n print(f\"-- SFTs amplitudes: {H1_SFTs.shape}\")\n H1_ts = f[ID_key]['H1']['timestamps_GPS']\n print(f\"-- timestamps: {H1_ts.shape} \\n\")\n\n # Retrieve the frequency data\n freq_data = np.array(f[ID_key]['frequency_Hz'])\n print(f\"- Frequency data: {freq_data.shape} \\n\")","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:47:54.952638Z","iopub.execute_input":"2022-10-26T21:47:54.953808Z","iopub.status.idle":"2022-10-26T21:47:54.990272Z","shell.execute_reply.started":"2022-10-26T21:47:54.953754Z","shell.execute_reply":"2022-10-26T21:47:54.989359Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load in the meta data ⏳","metadata":{}},{"cell_type":"code","source":"train_labels = pd.read_csv(DATA_PATH/'train_labels.csv')\ntrain_labels.head()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:47:54.99138Z","iopub.execute_input":"2022-10-26T21:47:54.991722Z","iopub.status.idle":"2022-10-26T21:47:55.012873Z","shell.execute_reply.started":"2022-10-26T21:47:54.99169Z","shell.execute_reply":"2022-10-26T21:47:55.011543Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Plot Train Test split ###\n\nplt.figure(figsize=(6,4))\nsns.barplot(['Train', 'Test'], [len(os.listdir(TRAIN_PATH)), len(os.listdir(TEST_PATH))]);\nplt.title(f'Train test split', fontsize=12)\nplt.ylabel('Count', fontsize=12)\nplt.xlabel('Category', fontsize=12)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:47:55.01449Z","iopub.execute_input":"2022-10-26T21:47:55.015252Z","iopub.status.idle":"2022-10-26T21:47:55.237571Z","shell.execute_reply.started":"2022-10-26T21:47:55.015209Z","shell.execute_reply":"2022-10-26T21:47:55.236607Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":" This is a very unusual train test split for a kaggle competetion. Normally we would see at an 80:20 split, but here it is 1:16. The competition hosts are encouraging us to generate more data ourselves ([see discussion](https://www.kaggle.com/competitions/g2net-detecting-continuous-gravitational-waves/discussion/347052))","metadata":{}},{"cell_type":"markdown","source":"## Labels\n Each data sample contains either **real or simulated noise** and possibly a **simulated continuous gravitational-wave signal** (CW). The task is to identify when a signal is present in the data (target=1)\n\n The `target_labels.csv` is a file containing the target labels; 1 if the data contains the presence of a gravitational wave, 0 otherwise. (Please note the presence of a small number of files labeled -1. Physicists are currently unable to determine the status of these files.)","metadata":{}},{"cell_type":"code","source":"labels_df = pd.read_csv(DATA_PATH/'train_labels.csv')\nlabels_df.head(5)","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:47:55.239261Z","iopub.execute_input":"2022-10-26T21:47:55.239881Z","iopub.status.idle":"2022-10-26T21:47:55.257Z","shell.execute_reply.started":"2022-10-26T21:47:55.239843Z","shell.execute_reply":"2022-10-26T21:47:55.255553Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"### Plot the distribution of labels ###\nlabel_count = labels_df['target'].value_counts()\nplt.figure(figsize=(10,8))\nsns.barplot(label_count.index, label_count.values, alpha=0.7)\nplt.title(f'Frequency of labels in training data', fontsize=12)\nplt.ylabel('Count', fontsize=12)\nplt.xlabel('label', fontsize=12)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:47:55.259713Z","iopub.execute_input":"2022-10-26T21:47:55.261429Z","iopub.status.idle":"2022-10-26T21:47:55.471114Z","shell.execute_reply.started":"2022-10-26T21:47:55.261357Z","shell.execute_reply":"2022-10-26T21:47:55.469695Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Sample submission:\n`sample_submission.csv` - a sample submission file in the correct format; your task is, for each file in the test/ folder, to predict the probability between [0, 1] of it it containing a continuous gravitational wave signal.","metadata":{"execution":{"iopub.status.busy":"2022-10-04T17:02:30.474744Z","iopub.execute_input":"2022-10-04T17:02:30.475143Z","iopub.status.idle":"2022-10-04T17:02:30.482773Z","shell.execute_reply.started":"2022-10-04T17:02:30.475112Z","shell.execute_reply":"2022-10-04T17:02:30.481314Z"}}},{"cell_type":"code","source":"pd.read_csv(DATA_PATH/'sample_submission.csv').head()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:47:55.473081Z","iopub.execute_input":"2022-10-26T21:47:55.473637Z","iopub.status.idle":"2022-10-26T21:47:55.502864Z","shell.execute_reply.started":"2022-10-26T21:47:55.473587Z","shell.execute_reply":"2022-10-26T21:47:55.501351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Plotting Spectograms 📊","metadata":{}},{"cell_type":"code","source":"### Helper functions to extract data and plot spectograms ###\n\ndef extract_data_from_hdf5(path, labels):\n '''\n Extracts data from hdf5 file and puts it into a dict. It also adds the label\n '''\n \n data = {}\n \n with h5py.File(path, \"r\") as f:\n\n ID_key = list(f.keys())[0]\n\n # Retrieve the frequency data\n data['freq'] = np.array(f[ID_key]['frequency_Hz'])\n\n # Retrieve the Livingston decector data\n data['L1_SFTs_amplitudes'] = np.array(f[ID_key]['L1']['SFTs'])\n data['L1_ts'] = np.array(f[ID_key]['L1']['timestamps_GPS'])\n\n # Retrieve the Livingston decector data\n data['H1_SFTs_amplitudes'] = np.array(f[ID_key]['H1']['SFTs'])\n data['H1_ts'] = np.array(f[ID_key]['H1']['timestamps_GPS'])\n \n # Get label from training labels if in training set\n data['label'] = labels.loc[labels.id==ID_key].target.item()\n \n return data\n \ndef plot_spectograms(data):\n '''\n Shows the real and imaginary amplitudes of the SFTs as spectograms for both detectors\n '''\n \n fig, ax = plt.subplots(2, 2, figsize=(16, 10))\n fig.suptitle(f\"Label {data['label']}\")\n\n for ind, detector in enumerate(['L1', 'H1']):\n ax[ind][0].set(xlabel=\"Timestamps [GPS]\",\n ylabel=\"Frequency [Hz]\",\n title=f\"{detector} - Real part\")\n ax[ind][1].set(xlabel=\"Timestamps [GPS]\",\n ylabel=\"Frequency [Hz]\",\n title=f\"{detector} - Imaginary part\")\n \n \n c0 = ax[ind][0].pcolormesh(data[f\"{detector}_ts\"], data['freq'],\n data[f\"{detector}_SFTs_amplitudes\"].real)\n c1 = ax[ind][1].pcolormesh(data[f\"{detector}_ts\"], data['freq'],\n data[f\"{detector}_SFTs_amplitudes\"].imag)\n \n fig.colorbar(c0, ax=ax[ind][0])\n fig.colorbar(c1, ax=ax[ind][1])\n \n plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:47:55.505115Z","iopub.execute_input":"2022-10-26T21:47:55.50617Z","iopub.status.idle":"2022-10-26T21:47:55.52446Z","shell.execute_reply.started":"2022-10-26T21:47:55.506114Z","shell.execute_reply":"2022-10-26T21:47:55.522886Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Lets plot some spectograms. One with a simulated signal and one without!","metadata":{}},{"cell_type":"code","source":"sns.reset_orig() # Reset seaborn theme (otherwise these plots come out red 😅)\ndata = extract_data_from_hdf5(train_example_with_signal_path, labels_df)\nplot_spectograms(data)","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:47:55.526034Z","iopub.execute_input":"2022-10-26T21:47:55.52646Z","iopub.status.idle":"2022-10-26T21:48:00.754122Z","shell.execute_reply.started":"2022-10-26T21:47:55.526424Z","shell.execute_reply":"2022-10-26T21:48:00.752351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data = extract_data_from_hdf5(train_example_without_signal_path, labels_df)\nplot_spectograms(data)","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:48:00.755731Z","iopub.execute_input":"2022-10-26T21:48:00.756608Z","iopub.status.idle":"2022-10-26T21:48:05.818948Z","shell.execute_reply.started":"2022-10-26T21:48:00.756565Z","shell.execute_reply":"2022-10-26T21:48:05.817593Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Very hard to notice any difference at all 🤔","metadata":{}},{"cell_type":"markdown","source":"# ","metadata":{}},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# What do the real and imaginary parts of a Short-time Fourier Transforms mean? 👻","metadata":{}},{"cell_type":"markdown","source":"## Plotting Spectograms\n\n## Basic spectrogram image classification \n COPIED FROM 🐢 JUN KODA\n\n","metadata":{}},{"cell_type":"markdown","source":"Power spectrum is averaged over every N = 32 timesteps:\n\n$$ P(\\omega) = \\frac{1}{N} \\sum_n^N \\left|a(\\omega, t_n) \\right|^2 $$\n\nImage size becomes: 2 channels x 360 frequencies x 128 timesteps.\n\nThe run time is dominated by data loading but I am not optimizing for that.\n\nVersion 4: Fix output target value to [0, 1] with `y_pred.sigmoid()`.","metadata":{}},{"cell_type":"code","source":"def power_spectrogram(h1_sft):\n img = np.empty((1,360, 128), dtype=np.float32)\n a=h1_sft[:, :4096] * 1e23\n p = a.real**2 + a.imag**2 # power\n p /= np.mean(p) # normalize\n p = np.mean(p.reshape(360, 128, 32), axis=2)\n img[0] = p\n fig, (ax1, ax2) = plt.subplots(1, 2,figsize=(5, 5))\n fig.suptitle('power density - and zoom of the zoom')\n ax1.imshow(img[0])\n ax2.imshow(img[0, 300:360])\n fig.tight_layout()\ndef power_spectrogram2(h1_sft):\n img = np.empty((1,360, 128), dtype=np.float32)\n a=h1_sft[:, :4096] * 1e23\n p = 2*(a.real**2 + a.imag**2)/(1800 * 1e24** 2) # power\n p /= np.mean(p) # normalize\n p = np.mean(p.reshape(360, 128, 32), axis=2)\n img[0] = p\n fig, (ax1, ax2) = plt.subplots(1, 2,figsize=(5, 5))\n fig.suptitle('power density - and zoom of the zoom')\n ax1.imshow(img[0])\n ax2.imshow(img[0, 300:360])\n fig.tight_layout()\n","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:48:05.82083Z","iopub.execute_input":"2022-10-26T21:48:05.821305Z","iopub.status.idle":"2022-10-26T21:48:05.83526Z","shell.execute_reply.started":"2022-10-26T21:48:05.821255Z","shell.execute_reply":"2022-10-26T21:48:05.833816Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"> **Variant 1: p = a.real\\*\\*2 + a.imag\\*\\*2 # power**","metadata":{}},{"cell_type":"code","source":"power_spectrogram(h1_sft)","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:48:05.837279Z","iopub.execute_input":"2022-10-26T21:48:05.837723Z","iopub.status.idle":"2022-10-26T21:48:06.16814Z","shell.execute_reply.started":"2022-10-26T21:48:05.837684Z","shell.execute_reply":"2022-10-26T21:48:06.166765Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":" > **Variant 2: 2\\*(a.real\\*\\*2 + a.imag\\*\\*2)/(1800 \\* 1e24\\*\\* 2) # normalized power**","metadata":{}},{"cell_type":"code","source":"power_spectrogram2(h1_sft)","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:48:06.173385Z","iopub.execute_input":"2022-10-26T21:48:06.173803Z","iopub.status.idle":"2022-10-26T21:48:06.563517Z","shell.execute_reply.started":"2022-10-26T21:48:06.173762Z","shell.execute_reply":"2022-10-26T21:48:06.56226Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## power 3D view","metadata":{}},{"cell_type":"code","source":"from mpl_toolkits.mplot3d import Axes3D\n#X = np.linspace(-20,20,100)\nX = np.linspace(0,360,num=360)\nY = np.linspace(0, 128, num=128)\nX, Y = np.meshgrid(X,Y)\nprint(X.shape,Y.shape)\na=h1_sft[:, :4096] * 1e24\np = np.absolute(a)+np.angle(a) # power\np /= np.mean(p) # normalize\np = np.mean(p.reshape(360, 128, 32), axis=2)\np.shape\n#Z = 4*xx**2 + yy**2\nZ=p.T\nfig = plt.figure()\nax = fig.add_subplot(111, projection='3d')\nax.plot_surface(X, Y, Z, cmap=\"hsv\", linewidth=0, antialiased=False, alpha=0.3)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:48:06.565442Z","iopub.execute_input":"2022-10-26T21:48:06.565816Z","iopub.status.idle":"2022-10-26T21:48:07.216011Z","shell.execute_reply.started":"2022-10-26T21:48:06.56578Z","shell.execute_reply":"2022-10-26T21:48:07.215076Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure(figsize=(8,8))\nax = fig.add_subplot(111, projection='3d')\nax.plot_surface(X, Y, Z, cmap=\"plasma\", linewidth=0, antialiased=False, alpha=0.2)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:50:28.752797Z","iopub.execute_input":"2022-10-26T21:50:28.753192Z","iopub.status.idle":"2022-10-26T21:50:29.567553Z","shell.execute_reply.started":"2022-10-26T21:50:28.753139Z","shell.execute_reply":"2022-10-26T21:50:29.566495Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Alternatives partial compression\nTo compress or to augument more?","metadata":{}},{"cell_type":"code","source":"img = np.empty((3,360, 128), dtype=np.float32)\na=h1_sft[:, :4096] * 1e23\np = 2*(a.real**2 + a.imag**2)/(1800 * 1e24** 2) # power\np /= np.mean(p) # normalize\np=p.reshape(360, 128, 32)\n#split the compression\np1=p[:,:,0:16]\np2=p[:,:,16:32]\np = np.mean(p, axis=2)\np1 = np.mean(p1, axis=2)\np2 = np.mean(p2, axis=2)\nimg[0] = p1\nimg[1] = p2\nimg[2] = p\nfig, (ax1, ax2,ax3) = plt.subplots(3, 1,figsize=(8, 8))\nplt.rcParams['image.cmap'] = 'BuPu'\nfig.suptitle('power density - and zoom of the zoom')\nax1.set_title(\"slice 0:16\")\nax2.set_title(\"slice 16:32\")\nax3.set_title(\"compressed 0:32\")\nax1.imshow(img[0, 300:360])\nax2.imshow(img[1, 300:360])\nax3.imshow(img[2, 300:360])\nfig.tight_layout()\n\n \n ","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:48:07.217595Z","iopub.execute_input":"2022-10-26T21:48:07.218473Z","iopub.status.idle":"2022-10-26T21:48:07.911232Z","shell.execute_reply.started":"2022-10-26T21:48:07.218432Z","shell.execute_reply":"2022-10-26T21:48:07.909887Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Timestamp analysis ⏱\n\nSince the continuous gravitational wave is simulated I am not sure how the length of it is determined (it seems to be hard coded in the [generating signals nb](https://github.com/PyFstat/PyFstat/blob/ec86602bb2f93238492a7242ad90995f6654eab7/examples/tutorials/1_generating_signals.ipynb)). As a result, I assume that the timestamp data is relatively meaningless... for completness I will include a few plots :)","metadata":{}},{"cell_type":"code","source":"### Extract timestamp data from training data ###\n\nH1_timestamps, L1_timestamps, start_diff, labels, freq = ([] for i in range(5))\n\nfor p in tqdm(os.listdir(TRAIN_PATH), total=len(os.listdir(TRAIN_PATH))):\n id_ = p.split('.')[0]\n labels.append(labels_df.loc[labels_df.id==id_].target.item())\n data = extract_data_from_hdf5(DATA_PATH/'train'/p, labels_df)\n L1_timestamps.append(data['L1_ts'])\n H1_timestamps.append(data['H1_ts'])\n start_diff.append(data['L1_ts'][0] - data['H1_ts'][0])\n freq.append(data['freq'])","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:48:07.913188Z","iopub.execute_input":"2022-10-26T21:48:07.913606Z","iopub.status.idle":"2022-10-26T21:50:26.17159Z","shell.execute_reply.started":"2022-10-26T21:48:07.913565Z","shell.execute_reply":"2022-10-26T21:50:26.169221Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = pd.DataFrame({'label':labels, 'L1_timestamp_length':[len(i) for i in L1_timestamps], 'H1_timestamp_length':[len(i) for i in H1_timestamps], 'Differnce in start time between detectors':start_diff})\ndf = df[df.label!=-1]","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:50:26.175845Z","iopub.execute_input":"2022-10-26T21:50:26.176343Z","iopub.status.idle":"2022-10-26T21:50:26.190544Z","shell.execute_reply.started":"2022-10-26T21:50:26.176294Z","shell.execute_reply":"2022-10-26T21:50:26.18919Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.set_theme()\n\nfig, ax = plt.subplots(1,3, figsize=(24,8))\nfig.suptitle(f\"In the plots the distribution of timestamps for both classes are shown; 1 indicates a simulated CW present and 0 not present\", fontsize=16)\nsns.histplot(\n df, x=\"L1_timestamp_length\", hue=\"label\",\n stat=\"density\", common_norm=False, bins=20, ax=ax[0], kde=True).set_title('Length of measurement for Livingston detector', fontsize=16);\n\nsns.histplot(\n df, x=\"H1_timestamp_length\", hue=\"label\",\n stat=\"density\", common_norm=False, bins=20, ax=ax[1], kde=True).set_title('Length of measurement for Hanford detector', fontsize=16);\n\nsns.histplot(\n df, x=\"Differnce in start time between detectors\", hue=\"label\",\n stat=\"density\", common_norm=False, bins=20, ax=ax[2], kde=True).set_title('Difference in starting timestamp between detectors', fontsize=16);","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:50:26.192398Z","iopub.execute_input":"2022-10-26T21:50:26.192763Z","iopub.status.idle":"2022-10-26T21:50:28.287025Z","shell.execute_reply.started":"2022-10-26T21:50:26.192731Z","shell.execute_reply":"2022-10-26T21:50:28.285701Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Frequency analysis 📳\n\nDo the detectors measure the same range of frequencies each time? No. The range varies but the minimum freq measured is ~50 Hz and the maximum is ~498 Hz","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(14,6))\nsns.histplot(x=list(np.hstack(freq)), binwidth=20)\nplt.title('Histogram of the range of Frequencies detected');\nplt.xlabel('Frequency Hz')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:50:28.28884Z","iopub.execute_input":"2022-10-26T21:50:28.289373Z","iopub.status.idle":"2022-10-26T21:50:28.714128Z","shell.execute_reply.started":"2022-10-26T21:50:28.289332Z","shell.execute_reply":"2022-10-26T21:50:28.712733Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Generating more data ⚙️","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Types of Noise\n- Non-stationary noise\n- Gaps\n- Narrow instrumental artifacts\n- Multiple detectors\n\nTypes of signal generation","metadata":{}},{"cell_type":"markdown","source":"# Baseline 📈","metadata":{}},{"cell_type":"markdown","source":"# Sample Submission","metadata":{}},{"cell_type":"code","source":"samp_sub = pd.read_csv(DATA_PATH/'sample_submission.csv')\nsamp_sub['target'] = 0.55\nsamp_sub.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2022-10-26T21:50:28.715611Z","iopub.execute_input":"2022-10-26T21:50:28.715939Z","iopub.status.idle":"2022-10-26T21:50:28.749561Z","shell.execute_reply.started":"2022-10-26T21:50:28.715908Z","shell.execute_reply":"2022-10-26T21:50:28.748442Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# WIP :)","metadata":{}},{"cell_type":"markdown","source":"#### Bibliography\n* https://kevinsprojects.wordpress.com/2014/12/13/short-time-fourier-transform-using-python-and-numpy/\n* https://www.ligo.org/science/GW-Continuous.php\n* https://link.springer.com/article/10.12942/lrr-2011-5\n* https://cplberry.com/2015/01/10/1408-0740/\n* https://www.aei.mpg.de/581560/a-milestone-in-the-race-to-find-the-first-continuous-gravitational-wave-signal\n* https://www.einstein-online.info/en/spotlight/continuousGW/","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}