?About the author: A Matlab simulation developer who loves scientific research. He cultivates his mind and improves his technology simultaneously.
For code acquisition, paper reproduction and scientific research simulation cooperation, please send a private message.
Personal homepage: Matlab Research Studio
Personal credo: Investigate things to gain knowledge.
For more complete Matlab code and simulation customization content, click
Intelligent optimization algorithm Neural network prediction Radar communication Wireless sensor Power system
Signal processing Image processing Path planning Cellular automaton Drone
Content introduction
Gravitational waves are predicted by Einstein’s general theory of relativity. They are fluctuations in space and time caused by extreme celestial motions. In 2015, gravitational waves were directly detected for the first time, a discovery considered a major breakthrough in physics. As gravitational wave detection technology continues to advance, we are now able to detect gravitational wave signals with unprecedented precision and frequency. However, the analysis and reproducibility of gravitational wave data is an important issue in gravitational wave research.
The repeatability of gravitational wave data means that different research teams can get similar results at different times and places. This characteristic is crucial for verifying and confirming the scientific significance of gravitational wave detection results. In gravitational wave research, data reproducibility can be ensured in a variety of ways.
First, standardization of data collection and processing is key to ensuring the repeatability of gravitational wave data. The operation and data collection of gravitational wave detectors require strict compliance with a series of standards and specifications to ensure the accuracy and comparability of data. At the same time, data processing and analysis processes also need to be standardized to ensure that different teams can obtain similar results when processing data.
Secondly, open data and open source code are important means to ensure the reproducibility of gravitational wave data. Data and code from gravitational wave research should be as open as possible so that other research teams can verify and replicate previous findings. Open data and open source code can also promote cooperation and exchanges in the field of gravitational wave research, helping to promote the development of the entire field.
In addition, the transparency and reproducibility of data analysis are also important factors in ensuring the repeatability of gravitational wave data. The research team should clearly document the process and methods of data analysis and provide sufficient information in the research results so that other teams can repeat the experiment following the same method. At the same time, scientific journals and academic publishing organizations should also encourage research teams to provide sufficient data and method details when publishing research results to promote data reproducibility.
Overall, the repeatability of gravitational wave data is an important issue in gravitational wave research. By standardizing data collection and processing procedures, opening data and open source code, and improving the transparency and reproducibility of data analysis, we can ensure the repeatability of gravitational wave data, thereby verifying and confirming the results of gravitational wave detection and promoting gravitational wave detection. Developments in the field of wave research. It is hoped that future gravitational wave research can make more progress in data reproducibility and reveal the mysteries of the universe for us.
Part of the code
%% Gravitational Wave Data Explorer</code><code>?</code><code>?</code><code>%% Introduction</code><code>% *Public Data:* Many public databases have been created for the purposes of </code><code>% making data freely accessible to the scientific community. A best practice is </code><code>% to assign a unique identifier to a dataset, so that it is discoverable. A common </code><code>% form of a unique identifier is a <https://en.wikipedia.org/wiki/Digital_object_identifier </code><code>% Digital Object Identifier or DOI> which points to the data. </code><code>% </code><code>% *Access Public Data:* To access and process public data, you can use several </code><code>% routes. </code> <code>?</code><code>% * Download data files to your local machine and work with them in MATLAB. </code><code>% * Access data directly via an API. MATLAB's <https://www .mathworks.com/help/matlab/ref/webread.html?searchHighlight=webread & amp;s_tid=srchtitle_webread_1 </code><code>% |webread|> function reads the RESTful API used by many portals.</code> <code>% * If the portal offers only Python bindings, <https://www.mathworks.com/help/matlab/call-python-libraries.html </code><code>% call Python from MATLAB>.</code><code>?</code><code>% *Data formats:* MATLAB supports a wide range of data formats</code><code>?</code><code>% * There are a wide range of scientific data formats that can be <https://www.mathworks.com/help/matlab/scientific-data.html </code><code>% natively read in MATLAB>. They include NetCDF and HDF5 as well as more specialized </code><code>% data formats. </code><code>% * In addition, MATLAB contains <https://www.mathworks.com/help/matlab/ref/fitsread.html </code> <code>% built-in functions> to read data in specific data formats often used in physics </code><code>% workflows.</code><code>% * Sometimes data import functions may be <https:// www.mathworks.com/matlabcentral/fileexchange/?category[]=support/sciences1689.support/physics1260 & amp;q=data + read </code><code>% written by the community>, and published on the MATLAB < https://www.mathworks.com/matlabcentral/fileexchange/ </code><code>% File Exchange> - a portal for community contributions in MATLAB. All community </code><code>% contributions are covered by open source licenses, which means they can be re-used, </code><code>% modified or added to. Exact terms and conditions depend on the licenses used </code><code>% by the author.</code><code>% * An example is the code used in this tutorial for strain data exploration </code><code>% and post-processing which follows this <https://de.mathworks.com/matlabcentral/fileexchange/108859-gravitationalwavedataexplorer ?tab=example & focused= </code><code>% community contribution>. In compliance with the Open Source License used for </code><code>% this code, the terms of the license are reprinted here. </code><code>?</code><code>% Copyright (c) 2022, Duncan Carlsmith All rights reserved. Redistribution and </code><code>% use in source and binary forms, with or without modification, are permitted </code><code>% provided that the following conditions are met: </code><code>% </code><code>% |* Redistributions of source code must retain the above copyright notice, </code><code>% this list of conditions and the following disclaimer.| </code><code>% </code><code>% |* Redistributions in binary form must reproduce the above copyright notice, </code><code>% this list of conditions and the following disclaimer in the documentation and/or </code><code>% other materials provided with the distribution| </code><code>% </code><code>% |* Neither the name of University of Wisconsin Madison nor the names of its </code><code>% contributors may be used to endorse or promote products derived from this software </code><code>% without specific prior written permission.| </code><code>% code><code>% </code><code>% |THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" </code><code>% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO , THE IMPLIED </code><code>% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. </code><code>% IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, </code> <code>% INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, </code><code>% BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, </code><code>% DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY </code><code>% OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE </code><code>% OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED </code><code>% OF THE POSSIBILITY OF SUCH DAMAGE.|</code><code>% </code><code>% </code> <code>% </code><code>% In this example, we will access gravitational wave data from the <https://www.gw-openscience.org/ </code><code>% GWOSC> Open Science database . </code><code>?</code><code>%% Access the data</code><code>% We will use the |webread| function to access the metadata using the RESTful </code><code>% API and use |h5read| function to read data directly from the url.</code><code>% Reading in archive of datafiles using API</code><code>% Clear workspace and specify URL endpoints</code> <code>?</code><code>baseURL = "https://www.gw-openscience.org/";</code><code>archive = webread(baseURL + "archive/all/json"); </code><code>?</code><code>% Read the archive of event files and sessions</code><code>archive.events;</code><code> </code><code>% Get event names</code><code>eventNames = string(fieldnames(archive.events));</code><code>writetable(table(eventNames),"eventNames.xlsx","WriteVariableNames",false)</code> code><code>?</code><code>% Select a given event</code><code>selectedEvent = eventNames(177);</code><code>eventName = extractBefore(selectedEvent,"_v"); </code><code>version = "v" + extract(selectedEvent,strlength(selectedEvent));</code><code>?</code><code>% Accessing metadata of selected event via API calls using |webread |</code><code>% Query the selected event</code><code>?</code><code>metadataEvent = webread(baseURL + "eventapi/json/event/" + erase(eventName,'x' ) + "/" + version);</code><code>?</code><code>% Show contents of the json</code><code>time_Event=metadataEvent.events.(selectedEvent).GPS;</code><code>?</code><code>% Show contents of the strain field</code><code>strainTable = struct2table(metadataEvent.events.(selectedEvent).strain);</code><code> ?</code><code>% Convert text to string:</code><code>strainTable.detector = string(strainTable.detector);</code><code>strainTable.format = categorical(strainTable.format); </code><code>strainTable.url = string(strainTable.url);</code><code>head(strainTable)</code><code>?</code><code>% Filter the table for only hdf5 files</code><code>hdf5Table = strainTable(strainTable.format == "hdf5",:);</code><code>% Use |h5read| within |webread| to read HDF5 data directly from the URL </code><code>% Select data to call (32 secs AND 4kHz sampling rate AND H1 or L1 detectors) </code><code>SF=4096;</code><code>available_detector_names=unique(hdf5Table.detector) ;</code><code>selectedDataset = hdf5Table(hdf5Table.duration == 32 & amp; ...</code><code> hdf5Table.sampling_rate == SF & amp; hdf5Table.detector == available_detector_names(1), :);</code><code>?</code><code>% The HDF5 file can be downloaded by entering the URL (in the last column of </code><code>% selectedDataset table) in the browser. However, we will not dowload the data, </code><code>% instead we will use |h5read| inside the |webread| function to read in data directly </code><code>% from the URL.</code><code>options = weboptions("ContentReader",@(x) h5read(x,"/strain/Strain"));</code><code>data = webread(selectedDataset.url,options);</code><code>?</code><code>% Plotting the raw (unfiltered) data</code><code>% Data just read is a time series data corresponding to unfiltered strain registered </code><code>% at one of the LIGO detectors. We normalize strain data and scale it to sensible </code><code>% values rather than 0.001 fm in 1 km.</code><code>figure</code><code>strain=( data-mean(data))*10^(19);</code><code>fs=selectedDataset.sampling_rate;</code><code>plotData(strain,fs,selectedDataset.detector)</code><code>saveas(gcf,"rawDataFigure.png")</code><code>?</code><code>%% Analyze the data</code><code>% Find noise peaks in the power spectrum using peak finder | findpeaks|</code><code>[p,~]=pspectrum(strain,fs);</code><code>q=10*log10(p);</code><code>?</code> <code>% Require peak height of 10 and sort in descending order of importance/height. </code><code>% We catch the peak values |pks|, locations |locs|, widths |w|, and prominences </code><code>% code><code>% |prom|. </code><code>% </code><code>% Call |findpeaks| again without collecting results to make a plot.</code><code>figure</code><code>findpeaks(q,MinPeakProminence=5, SortStr='descend');</code><code>title(strcat(selectedDataset.detector,' PSD with noise peaks identified by FindPeaks'));</code> <code>xlabel('frequency (Hz)');</code><code>ylabel('Power Spectrum of Raw Strain (dB)');</code><code>saveas(gcf,"noisePeaksSpectrumFigure.png" )</code><code>?</code><code>%% Apply notch and bandpass filters to remove noise and plot the power spectrum of the filtered data</code><code>% We will next apply notch and bandpass filter with lower limit 35Hz and upper </code><code>% limit 350Hz (corresponding to the relevant frequency range for detecting first </code><code>% gravitational wave by the LIGO detectors) as in published analysis [1]. A notch </code><code>% filter suppresses frequencies within a narrow range and we apply it to remove </code><code>% various peaks previously found, since most pronounced peaks most likely correspond </code><code> % to noise that we would like to remove. See function |filterSignal| at the end </code><code>% of this script. </code><code>strainfilt=filterSignal(strain,fs);</code> <code>figure</code><code>pspectrum(strainfilt,fs);</code><code>title(selectedDataset.detector + ' filtered strain PSD');</code><code>subtitle(eventName)</code><code>saveas(gcf,"noiseFilteredSpectrumFigure.png")</code><code>?</code><code>%% Plot the filtered strain data in a 0.3s time-window around the event</code><code> code><code>% Standard deviation of the strain signal.</code><code>strainsig=std(strainfilt);</code><code>?</code><code>% Create a 0.3s time window around published time near the center of the 32 </code><code>% s long signal,</code><code>tstart=time_Event-selectedDataset.GPSstart-0.15;</code><code>tend=tstart + 0.3; </code><code>time=(1:length(strain))/fs;</code><code>time=time';</code><code>?</code><code>% Make initial mask to drop endpoints and glitches from data.</code><code>mask=time<tend & time>tstart & strainfilt<3*strainsig & strainfilt>-3*strainsig;</code><code> </code><code>% Plot the raw and filtered subsamples.</code><code>figure</code><code>tiledlayout(2,1);</code><code>nexttile</code><code>plot(time(mask),strain(mask))</code><code>ylabel('Raw Strain (\times 10^{19})')</code><code>title("Strain timeseries data centered at the event")</code><code>subtitle(eventName)</code><code>legend(selectedDataset.detector,Location='best')</code><code>grid on</code><code>nexttile</code><code>plot(time(mask),strainfilt(mask))</code><code>xlabel('Time(s)')</code><code>ylabel(' Filtered Strain (\times 10^{19})')</code><code>legend(selectedDataset.detector,Location='best')</code><code>grid on</code><code>saveas( gcf,"filteredStrainDataFigure.png")</code><code>?</code><code>%% Make spectrogram of chirp</code><code>% As the black holes inspiral and orbit faster and faster before merging in </code><code>% one black hole, the frequency of the gravitational radiation increases. A spectrogram </code><code>% exhibits the time variation of a power spectrum. At any given time, the spectrum </code> <code>% is derived from window around that time. We will obtain a time-frequency representation </code><code>% of the filtered strain data, showing the signal frequency increasing over time. </code><code>% </code><code>% Make a spectrogram using |pspectrum| with the specifier 'spectrogram' that </code><code>% shows the chirp of frequency. Specify to a time window width over which to </code><code>% calculate the spectrum, a width small compared to the duration of the signal, </code><code>% and a large fractional overlap of successive windows.</code><code>figure</code><code> pspectrum(strainfilt(mask),time(mask),'spectrogram',FrequencyLimits=[32,400], ...</code><code> TimeResolution=0.033,OverlapPercent=97,Leakage=0.5);</code><code>title("Time-frequency representation of the strain data")</code><code>subtitle(eventName)</code><code>saveas(gcf,"chirpSpectrogramFigure.png")</code><code> %text(time_Event-selectedDataset.GPSstart + 0.015,250,strcat(selectedDataset.detector,' chirp \rightarrow'),HorizontalAlignment='right');</code><code>?</code><code>% The chirp is most evident as the rising band between frequencies 50Hz to 300Hz.</code><code>?</code><code>%% Compute correlation between the two detectors</code><code>% Let's bring in the strain data from the same gravitational wave event, but </code><code>% from another interferometer, into data2 variable and filter out signal as we </code><code>% did above.</code><code>switch selectedDataset .detector</code><code> case 'H1'</code><code> selectedDataset2 = hdf5Table(hdf5Table.duration == 32 & hdf5Table.sampling_rate == SF & hdf5Table.detector=='L1' ,:);</code><code> if isempty(selectedDataset2.detector)</code><code> selectedDataset2 = hdf5Table(hdf5Table.duration == 32 & amp; hdf5Table.sampling_rate == SF & amp; hdf5Table.detector =='V1',:);</code><code> end</code><code> otherwise</code><code> selectedDataset2 = hdf5Table(hdf5Table.duration == 32 & amp; hdf5Table.sampling_rate == SF & amp; hdf5Table.detector=='H1',:);</code><code>end</code><code>?</code><code>data2= webread(selectedDataset2.url,options); </code><code>strain2=(data2-mean(data2))*10^(19);</code><code>strain2filt=filterSignal(strain2,fs);</code><code>?</code> code><code>% Now we calculate correlation between the two filtered signals at two LIGO </code><code>% interferometers separated by ~3,000km straight line distance from each other.</code><code>[acor,lag ] = xcorr(strainfilt(mask),strain2filt(mask));</code><code>?</code><code>% Find the maximum in the correlation and the corresponding time difference.</code><code> [~,I] = max(abs(acor));</code><code>lagDiff = lag(I);</code><code>timeDiff = lagDiff/fs;</code><code>?</code><code>% The time lag of $\sim 7$ms for the event GW150914 is within the range of the </code><code>% quoted value of $6.9^{ + 0.5}_{-0.4}$ ms [1].</code><code>% </code><code>% Detection of the signal separately by each of the two LIGO detectors (both </code><code>% of which compare favorably with the simulated strain from general relativity </code><code>% shown below), within a time window of $\pm10$ms was crucial to claim the first </code><code>% ever direct experimental observation of the gravitational wave. </code><code>figure</code><code>plot(lag/fs,acor)</code><code>title('Time Correlation of ' + selectedDataset.detector + ' and ' + selectedDataset2.detector + ' Signals' )</code><code>subtitle(eventName)</code><code>xlabel('Time (s)');</code><code>ylabel('Correlation (arb. units)')</code><code>grid on</code><code>saveas(gcf,"detectorCorrelationFigure.png")</code><code>%% Plot two signals at two detectors together</code><code>% For comparing two signals with each other, following LIGO collaboration result </code><code>% [1], we multiply one of the filtered strains with |sign(acor(I))| and add time </code><code> % lag to it.</code><code>?</code><code>plot(time(mask),strainfilt(mask),time(mask) + timeDiff,sign(acor(I))*strain2filt(mask ))</code><code>xlabel('Time(s)')</code><code>ylabel('Filtered Strains (\times 10^{19})')</code><code>title( selectedDataset.detector + ' and ' + selectedDataset2.detector + ' Signals')</code><code>subtitle(eventName)</code><code>xlim([tstart,tend])</code><code>legend (selectedDataset.detector,selectedDataset2.detector,Location='best');</code><code>?</code><code>?</code><code>%% References</code><code>% [1] B. P. Abbott, Benjamin _et al._ (LIGO Scientific Collaboration and Virgo </code><code>% Collaboration) "Observation of Gravitational Waves from a Binary Black Hole </code><code>% Merger", <https ://journals.aps.org/prl/abstract/10.1103/PhysRevLett.116.061102 </code><code>% _Phys. Rev. Lett._ *116,* 061102> (2016)</code><code>% </code><code>% [2] R. Abbott _et al_. (LIGO Scientific Collaboration, Virgo Collaboration </code><code>% and KAGRA Collaboration), "Open data from the third observing run of LIGO, Virgo, </code><code>% KAGRA and GEO", <https://doi.org/10.3847/1538-4365/acdc9f ApJS 267 29 (2023)> </code><code>% -- <https:/ /inspirehep.net/literature/2630216 INSPIRE> </code><code>% </code><code>% [3] R. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), </code><code>% "Open data from the first and second observing runs of Advanced LIGO and Advanced </code><code>% Virgo", <https://doi.org/10.1016/j.softx.2021.100658 SoftwareX 13 (2021) 100658> </code><code>% -- <https://inspirehep.net/literature/1773351 INSPIRE></code><code>?</code><code>%% Acknowledgements</code><code>%This research has made use of data or software obtained from the Gravitational Wave </code><code>% Open Science Center (gwosc.org), a service of the LIGO Scientific Collaboration, </code><code>% the Virgo Collaboration, and KAGRA. This material is based upon work supported by </code><code>% NSF's LIGO Laboratory which is a major facility fully funded by the National Science </code><code>% Foundation, as well as the Science and Technology Facilities Council (STFC) of the </code><code>% United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany </code><code>% for support of the Construction of Advanced LIGO and construction and operation of the </code><code>% GEO600 detector. Additional support for Advanced LIGO was provided by the Australian </code><code>% Research Council. Virgo is funded, through the European Gravitational Observatory (EGO), </code><code>% by the French Center National de Recherche Scientifique (CNRS), the Italian Istituto </code><code>% Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by </code><code>% institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan, Monaco, Poland, </code><code>% Portugal, Spain. KAGRA is supported by Ministry of Education, Culture, Sports , Science </code><code>% and Technology (MEXT), Japan Society for the Promotion of Science (JSPS) in Japan; </code><code>% National Research Foundation (NRF) and Ministry of Science and ICT ( MSIT) in Korea; </code><code>% Academia Sinica (AS) and National Science and Technology Council (NSTC) in Taiwan</code><code>?
Operation results
