基于Matlab演示引力波数据可重复的分析

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