Skip to content
Snippets Groups Projects
Commit ccb7139a authored by Stefano Covino's avatar Stefano Covino
Browse files

120325 commit

parent 6fadd8cb
Branches
Tags
No related merge requests found
Showing
with 912 additions and 330 deletions
......@@ -109,7 +109,7 @@
"provenance": []
},
"kernelspec": {
"display_name": "Julia 1.11.3",
"display_name": "Julia 1.11.4",
"language": "julia",
"name": "julia-1.11"
},
......@@ -117,7 +117,7 @@
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.11.3"
"version": "1.11.4"
}
},
"nbformat": 4,
......
%% Cell type:markdown id:91330533 tags:
**What is this?**
*This jupyter notebook is part of a collection of notebooks on various topics discussed during the Time Domain Astrophysics course delivered by Stefano Covino at the [Università dell'Insubria](https://www.uninsubria.eu/) in Como (Italy). Please direct questions and suggestions to [stefano.covino@inaf.it](mailto:stefano.covino@inaf.it).*
%% Cell type:markdown id:915ee876 tags:
**This is a textual notebook**
%% Cell type:markdown id:53194e25 tags:
![Time Domain Astrophysics](Pics/TimeDomainBanner.jpg)
%% Cell type:markdown id:a7b36f9a tags:
# Introduction
***
%% Cell type:markdown id:b336520e tags:
## Contacts
***
![Stefano](Pics/Stefano.png)
- Stefano Covino
- INAF / Brera Astronomical Observatory
- +39 02 72320475
- +39 3316748534 (if urgent…)
- Emails: [stefano.covino@inaf.it](mailto:stefano.covino@inaf.it) - [stefano.covino@uninsubria.it](mailto:stefano.covino@uninsubria.it)
- Web: https://sites.google.com/a/inaf.it/stefano-s-site/
![Banner](Pics/Banner.png)
%% Cell type:markdown id:9556e54d tags:
## Main Goal of the course: Have fun!
***
|![data](Pics/data.jpg)|![regression](Pics/regression.jpg)|
|----------------------|----------------------------------|
%% Cell type:markdown id:cffcb879 tags:
## Program (for 6 or 7 courses, roughly…)
***
1. Introduction to time series
2. Time (and spatial) variability in astrophysics
3. Fourier analysis and noise characterization
4. Case study: stellar variability
5. Case study: exo-planet transits
6. Case study: pulsars
7. Time-domain analysis and auto-regressive processe
8. Irregular sampling, Lomb-Scargle periodograms
9. Case studies: AGN variability
10. Advanced topics: non-parametric analysis
11. Matching filters
12. Case study: LIGO/Virgo gravitational wave signals
13. Data exploration
14. Case study: SETI data analysis
15. Big-data, machine learning and “intelligent” systems for time-series analysis
16. Case studies: spatial variability (CMB, large scale structure)
17. Final topics: forecasting
> In reality these are just topics that can be covered. We can stress different aspects depending on the interests of the *students*.
%% Cell type:markdown id:2d351400 tags:
## Time-Series are ubiquitous
***
- Anytime we have a measurement repetated multiple times we have a time-series.
![CO_2 vs Temp](Pics/CO2T.png)
![CO_2](Pics/CO2.png)
![Neptune](Pics/Neptune.png)
- As a matter of fact, a time-series does not need to have "time" as index!
![PAMELA](Pics/PAMELA.png)
![Satellite](Pics/satellite.png)
%% Cell type:markdown id:722fb437 tags:
## Temptative schedule (don’t trust it too much…)
***
1. 26/2 - Introduction
2. 27/2 - Statistics reminder - part I
3. 5/3 - Statistics reminder - part II
4. 6/3 - Spectral analysis - part I
5. 12/3 - Spectral analysis - part II
6. 13/3 - Science cases: Sunspots Number - X-ray Binaries
7. 19/3 - Irregularly sampled time series - part I
8. 20/3 - Irregularly sampled time series - part II
9. 26/3 - Science Cases - Variable Stars - AGN and blazars
10. 27/3 - Time domain analysis - part I
11. 2/4 - Time domain analysis - part II
12. 3/4 - Guest lecture - Spectral analysis in Cosmology
13. 9/4 - Guest lecture - X-ray pulsators
14. 10/4 - Time domain analysis - ARIMA models
15. 16/4 - Time domain analysis - Advanced tools
16. 30/4 - Wavelet analysis
17. 7/5 - Guest lecture - Exoplanets
18. 8/5 - Time of arrival analysis
19. 14/5 - Non-parametric methods
20. 15/5 - Gaussian processes
21. 21/5 - Science case: GRBs
22. 22/5 - Astrostatistics final considerations
%% Cell type:markdown id:f73fd0ea tags:
## How is the course managed?
***
### Frontal lectures
- These are the traditional university lectures.
- Although this increases the organizational complexity substantially, I am availbale to stream and record my lectures, if needed.
- There are contraindications. As a matter of fact, this is one of few cases where a remote access is not even close as effective as being in presence.
![Frontal Lectures](Pics/FrontalLectures.jpg)
### Real research life examples…
- Scientists working in the field will deliver "didactic lectures", allowing one to see most of ideas deveooped during the course applied in a real research environment.
![Paperino astronomo](Pics/Paperino.jpg)
### (Optional) papers to deepen our knowledge…
- Most of the topics discussd during the course can be investigated thoroughly and papers from astrophysical (mainly) literature are presented for particularly concerned readers.
![Author et al.](Pics/Papersetal.jpg)
### Question time
- The course is divided in several main sections. At the end of each of them, some time will be devoted to open discussions and questions.
![Questions](Pics/Questions.gif)
### Lectures from specialists in the field
- Together with regular lectures, a few specialists in the field, i.e. scientist carrying out researches by time-domain tools and techniques, are invited to describe their works.
![Nilus](Pics/Nilus.jpg)
### Language
- According to university guidelines, lectures will be delivered in English. Of course, a fair evaluation of the context might ask some flexibility.
![Language](Pics/language.jpg)
### Statistical framework
- During this course we are going to work in a Bayesian framework.
- Bayesian statistics is an approach to inferential statistics based on Bayes' theorem, where available knowledge about parameters in a statistical model is updated with the information in observed data. The background knowledge is expressed as a prior distribution and combined with observational data in the form of a likelihood function to determine the posterior distribution. The posterior can also be used for making predictions about future events.
- Nevertheless, we are not dogmatic and mentions or applications based on familiar "frequentist" approaches are preseneted and discussed, when we deem it opportune.
![Statistics](Pics/Bayesians.png)
### Programming languages
- Most of the examples we are going to analyze during the course are based on some sort of computer analysis.
- `Python` is *de-facto* the standard language in data science.
- Yet, while this language is definitely truly amazing, well designed and worth mastering, for the specific needs of scientific computing there are alternatives of growing popularity.
- We threfore provide examples mainly with `Julia`, and encourage the students to get some confidence with this programming language too.
- Notebooks are written by the [markdown language](https://www.markdownguide.org/basic-syntax/), a simple language integrating features of the HTML and latex languages.
|![python](Pics/python.png)|![julia](Pics/julia.png)|
|--------------------------|------------------------|
%% Cell type:markdown id:0a132f61-6a28-4c74-a63d-b1442f23269c tags:
- A remarkable introducti0n to the `julia` language for a scientist is available online [here](https://juliadatascience.io/).
%% Cell type:markdown id:9bf3b78e tags:
## Warning! The course is not only for astrophysicists!
- It is indeed part of the set of courses for future astrophysicists. Nevertheles, almost nothing we are going to discuss is truly only for astrophysics. In reality, several applications and ideas are taken from other fields, i.e. economics, social sciences, climatology, etc.
![Physicists and astrophysicists](Pics/astrophysics.jpg)
%% Cell type:markdown id:dedc01e8 tags:
## Final assessment
- The final examination is an oral one.
- *Students* must interact with the teacher in advance of the examination and a science case obtained by the modern literature will be selected.
- The *student* will be asked to properly describe the main formal aspects of the study and discuss critically the reliability and limits of the presented results.
%% Cell type:markdown id:6e9e2d08 tags:
## Gitlab repository
- Slides, notebooks, papers, etc. are available on [gitlab](https://www.ict.inaf.it/gitlab/stefano.covino/TimeDomainAstrophysics.git)
- Check the repository frequently since is (rather often) updated during the course.
|![gitlab](Pics/gitlab.jpg)|![course](Pics/gitlabcourse.png)|
|-------------------------|-------------------------------|
%% Cell type:markdown id:e96ac78c tags:
## Relaxing time(-series...)
![Relax](Pics/relaxing.png)
%% Cell type:markdown id:ffe137ba tags:
## Reference & Material
- The course is based on published scientific papers distributed by the teacher before any main topic is addressed.
- Science cases are based on actual scientific papers as well.
- Slides prepared by the teacher will also be distributed.
- A general introductory text to time series analysis as: [“Introduction to Time Series and Forecasting”, by P.J. Brockwell and R.A Davis](https://link.springer.com/book/10.1007/978-3-319-29854-2) might be useful. However, any other analogous text easily obtainable by the student will be fine as well.
- Two textbooks more strictly related to the topics discussed during the course mainly, but not only, for astrophysical applications are:
- [“Modern Statistical Methods for Astronomy”, by E.D. Feigelson and G.J. Babu](https://www.cambridge.org/core/books/modern-statistical-methods-for-astronomy/941AE392A553D68DD7B02491BB66DDEC)
- [“Statistics, data Mining and Machine Learning in Astronomy”, by Ivezić et al.](https://press.princeton.edu/books/hardcover/9780691198309/statistics-data-mining-and-machine-learning-in-astronomy)
%% Cell type:markdown id:dff2c952 tags:
## Further Material
Papers for examining more closely some of the discussed topics.
- [Voughan et al. (2013) - "Random Time Series in Astronomy"](https://royalsocietypublishing.org/doi/10.1098/rsta.2011.0549)
- [Storopoli et al. (2021) - "Julia Data Science"](https://juliadatascience.io/)
%% Cell type:markdown id:05e93b1d tags:
## Course Flow
<table>
<tr>
<td>Previous lecture</td>
<td>Next lecture</td>
</tr>
<tr>
<td><a href="../../Course.ipynb">Course Summary</a></td>
<td><a href="../Lecture%20-%20Statistics%20Reminder/Lecture-StatisticsReminder.ipynb">Statistics Reminder</a></td>
</tr>
</table>
%% Cell type:markdown id:591bd355 tags:
**Copyright**
This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the examples, unless obtained from other properly quoted sources, under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Stefano Covino, Time Domain Astrophysics - Lecture notes featuring computational examples, 2025*.
......
......@@ -273,6 +273,7 @@ md"""
Papers for examining more closely some of the discussed topics.
- [Voughan et al. (2013) - "Random Time Series in Astronomy"](https://royalsocietypublishing.org/doi/10.1098/rsta.2011.0549)
- [Storopoli et al. (2021) - "Julia Data Science"](https://juliadatascience.io/)
"""
# ╔═╡ c99cc57a-d012-4d11-a3e6-1e2ba35ce92e
......@@ -314,7 +315,7 @@ PlutoUI = "~0.7.61"
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised
julia_version = "1.11.3"
julia_version = "1.11.4"
manifest_format = "2.0"
project_hash = "6d1b77f27e79835fc27b2d7e99ab8fcaf37aa976"
......
......@@ -1450,7 +1450,7 @@
"provenance": []
},
"kernelspec": {
"display_name": "Julia 1.11.3",
"display_name": "Julia 1.11.4",
"language": "julia",
"name": "julia-1.11"
},
......@@ -1458,7 +1458,7 @@
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.11.3"
"version": "1.11.4"
}
},
"nbformat": 4,
......@@ -747,14 +747,23 @@ md"""
"""
# ╔═╡ 1112b278-802d-4611-8954-2795db178e88
pf = @bind pl_fr NumberField(1:10, default=2);
begin
pf = @bind pl_fr NumberField(1:10, default=2)
af = @bind am_fr NumberField(100:100:1000, default=300)
end;
# ╔═╡ 19a6c251-adc1-422d-a702-20ea9b971449
cm"Frequency for the test plot:"
cm"Frequency for the test signal:"
# ╔═╡ a8b99222-18a7-4b47-8f98-8961a708ec57
pf
# ╔═╡ 77b0596d-8c4a-4d3d-9bd8-54d3fde618c8
cm"Amplitude for the test signal:"
# ╔═╡ d542c639-1c93-4edd-9132-facd18fc89fb
af
# ╔═╡ b96525bf-90e5-47b3-b305-1db0ad675c2c
begin
dt = 0.03125 # seconds
......@@ -763,8 +772,8 @@ begin
times = range(start=0, stop=exposure-dt, step=dt) # seconds
long_times = range(start=0, stop=long_exposure-dt, step=dt) # seconds
signal = 300 .* sin.(pl_fr .* pi .* times ./ 0.5) .+ 1000 # counts/s
long_signal = 300 .* sin.(pl_fr .* pi .* long_times ./ 0.5) .+ 1000 # counts/s
signal = am_fr .* sin.(pl_fr .* pi .* times ./ 0.5) .+ 1000 # counts/s
long_signal = am_fr .* sin.(pl_fr .* pi .* long_times ./ 0.5) .+ 1000 # counts/s
noisy = [rand(Poisson(theta)) for theta in signal .* dt] # counts
long_noisy = [rand(Poisson(theta)) for theta in long_signal .* dt] # counts
......@@ -792,6 +801,9 @@ md"""
# ╔═╡ b7213b38-21ac-4fb5-ab53-694d3ac2fda9
pf
# ╔═╡ 78dc57ce-1cbd-4f9c-9741-54bf430b9186
af
# ╔═╡ 25fcf199-7da4-4632-a2fa-d0131f58d533
begin
# We use the DSP package
......@@ -1233,10 +1245,13 @@ This notebook is provided as [Open Educational Resource](https://en.wikipedia.or
# ╟─b2e2bace-7e50-41f0-87fb-780737ff6616
# ╟─1112b278-802d-4611-8954-2795db178e88
# ╠═19a6c251-adc1-422d-a702-20ea9b971449
# ╟─a8b99222-18a7-4b47-8f98-8961a708ec57
# ╠═a8b99222-18a7-4b47-8f98-8961a708ec57
# ╟─77b0596d-8c4a-4d3d-9bd8-54d3fde618c8
# ╟─d542c639-1c93-4edd-9132-facd18fc89fb
# ╟─b96525bf-90e5-47b3-b305-1db0ad675c2c
# ╟─ae8a4ebf-da8d-4778-ad9a-87e3aab72399
# ╟─b7213b38-21ac-4fb5-ab53-694d3ac2fda9
# ╟─78dc57ce-1cbd-4f9c-9741-54bf430b9186
# ╟─25fcf199-7da4-4632-a2fa-d0131f58d533
# ╟─3e27b7d8-9dec-4e54-9050-5ca9275d8499
# ╟─817993c0-a1ce-46d2-a213-ed6448d5159f
......
%% Cell type:markdown id:91330533 tags:
 
**What is this?**
 
 
*This jupyter notebook is part of a collection of notebooks on various topics discussed during the Time Domain Astrophysics course delivered by Stefano Covino at the [Università dell'Insubria](https://www.uninsubria.eu/) in Como (Italy). Please direct questions and suggestions to [stefano.covino@inaf.it](mailto:stefano.covino@inaf.it).*
 
%% Cell type:markdown id:915ee876 tags:
 
**This is a `Julia` notebook**
 
%% Cell type:code id:e29ceecb tags:
 
``` julia
import Pkg; Pkg.activate(".")
```
 
%% Output
 
 Activating project at `/mnt/chromeos/GoogleDrive/MyDrive/Lab/Teaching/Insubria/Docs_2024_25/Lectures/Lecture - Wavelet Analysis`
 
%% Cell type:code id:08b22357 tags:
 
``` julia
Pkg.instantiate()
```
 
%% Cell type:code id:4ff28a4e tags:
 
``` julia
#using AbstractFFTs
using CairoMakie
#using CSV
using ContinuousWavelets
#using DataFrames
#using Dates
using Distributions
using DSP
#using FFTW
#using FITSIO
#using Format
#using HTTP
#using HypothesisTests
#using LaTeXStrings
#using LinearAlgebra
#using LombScargle
#using Optim
#using PValue
using Random
#using Statistics
```
 
%% Cell type:markdown id:53194e25 tags:
 
![Time Domain Astrophysics](Pics/TimeDomainBanner.jpg)
 
%% Cell type:markdown id:a7b36f9a tags:
 
# Wavelet Analysis
 
***
 
## PSD Varying Time-Series
***
 
- The trigonometric basis functions used in the Fourier transform have an infinite extent and for this reason the Fourier transform may not be the best method to analyze non periodic time series data, such as the case of a localized event (e.g., a burst that decays over some time scale so that the PSD is also varying with time).
 
- We can (and actually do) evaluate the PSD for finite stretches of time series in order to detect its changes.
 
- This approach (called *spectrogram*, or *dynamical power spectra analysis*, or even *windowed Fourier transform*) suffers from degraded spectral resolution and is sensitive to the specific choice of time series segmentation length.
 
### Windowed Fourier transform
***
 
- The FT is performed on a sliding segment of length $T$ from a time series of time step $δt$ and total length $Nδt$, thus returning frequencies from $T^{−1}$ to $2δt^{−1}$ at each time step.
 
- The segments can be windowed with an arbitrary function such as a boxcar (no smoothing) or any other choice.
 
- The WFT represents an inaccurate and inefficient method of time–frequency localization, as it imposes a scale or “response interval” T into the analysis.
 
- The inaccuracy arises from the aliasing of high- and low-frequency components that do not fall within the frequency range of the window. The inefficiency comes from the $T/(2δt)$ frequencies, which must be analyzed at each time step.
 
 
%% Cell type:markdown id:5025348f tags:
 
## Wavelets
***
 
- A popular family of basis functions is called wavelets.
 
- By construction, wavelets are localized in both frequency and time domains.
 
- Individual wavelets are specified by a set of *wavelet filter coefficients*. Given a wavelet, a complete orthonormal set of basis functions can be constructed by scalings and translations.
 
- Different wavelet families trade the localization of a wavelet with its smoothness. Popular wavelets include “Mexican hat”, Haar and Daubechies wavelets.
 
- A wavelet is a “small wave,” where small derives from the fact that it is mostly limited to an interval of time.
 
- A (mother) wavelet $ψ$ has to satisfy two requirements: its integral must be zero and the integral of its square must be unity:
 
$$ \int_{-\infty}^\infty \Psi(t) dt = 0 \qquad \int_{-\infty}^\infty \Psi^2(t) dt = 1 $$
 
- It is easy to see that these requirements imply that $ψ(t)$ is essentially non-zero only over a limited range of t and that it has to extend both above and below zero.
 
![Wavelets](Pics/wavelets.png)
 
- Three examples of wavelets. From left to right: Haar wavelet, mexican hat and morlet wavelet.
 
%% Cell type:markdown id:e9058425 tags:
 
- A wavelet can be shifted by $τ$ and dilated by a scale parameter $σ$:
 
$$ \Psi_{\tau,\sigma}(t) = \frac{1}{\sqrt{\sigma}} \Psi\left(\frac{t-\tau}{\sigma}\right) $$
 
- Again $ψ(t)$ is essentially non-zero only over a limited range of $t$ and that it has to extend both above and below zero.
 
<br>
 
- The wavelet transform of a function $f(t)$ is computed by correlating $f(t)$ with the complex conjugate of $ψ_{τ,σ}(t)$ (wavelets can also be complex functions):
 
$$ W[f(\tau,\sigma)] = \int_{-\infty}^\infty f(t) \Psi_{\tau,\sigma}^*(t)dt$$
 
- The *discrete wavelet transform* (DWT) can be used to analyze the power spectrum of a time series as a function of time.
 
- A possible wavelet might be: $w(t|t_0, f_0, Q) = A \exp[i 2 \pi f_0 (t-t_0)] \exp[-f_0^2 (t-t_0)^2 / Q^2]$,
- where $t_0$ is the central time, $f_0$ is the central frequency, and the dimensionless parameter $Q$ is a model parameter which controls the width of the frequency window.
 
![Wavelet example](Pics/waveletex1.png)
 
- The Fourier transform of the previous wavelet is:
 
$$ W(t|t_0, f_0, Q) = \left( \frac{\pi}{f_0^2/Q^2} \right)^{1/2} \exp[i 2 \pi f_0 t_0] \exp[-\frac{\pi^2 Q^2 (f-f_0)^2}{Q f_0^2}]$$
 
- Technically speaking, this is not exactly a wavelet. It is indeed closely related to a true wavelet, the Morlet wavelet, through a simple scaling and offset.
- One might think to it as a “matched filters”.
 
- Given a signal $h(t)$, its wavelet transform is given by the convolution:
 
$$ H_w(t_0; f_0, Q) = \int_{-\infty}^\infty h(t)w(t|t_0, f_0, Q) $$
 
- By the convolution theorem, we can write the Fourier transform of $H_w$ as the pointwise product of the Fourier transforms of $h(t)$ (e.g. by a DFT) and of $w(t | t_0, f_0, Q)$.
 
- Then, the wavelet PSD, can be defined by:
 
$$ {\rm PSD}_w(f_0, t_0; Q) = |H_w(t_0; f_0, Q)|^2 $$
 
- As said, unlike the typical Fourier-transform PSD, the wavelet PSD allows detection of frequency information which is localized in time.
 
- This is indeed one the approaches (we are simplifying a lot here…) used in the, e.g., LIGO/Virgo projects to detect gravitational wave events.
 
- Because of the noise level in the GW measurements, rather than a standard wavelet they instead use functions which are tuned to the expected form of the signal (i.e., matched filters).
 
<br>
 
- Example with a localized Gaussian noise:
 
![Wavelet example](Pics/waveletex2.png)
 
- The upper panel shows the signal. The middle panel shows an example wavelet and the lower panel shows the power spectral density as a function of the frequency $f_0$ and the time $t_0$, for $Q = 1.0$.
 
<br>
 
- And now a different example:
 
![Wavelet example](Pics/waveletex3.png)
 
- The upper panel shows the input signal, which consists of a Gaussian spike in the presence of white (Gaussian) noise. The middle panel shows again an example wavelet. The lower panel shows the PSD as a function of the frequency $f_0$ and the time $t_0$, for $Q = 0.3$.
 
%% Cell type:markdown id:d4abc60f tags:
 
#### Exercize: wavelet analysis of a chirp signal
***
 
- Let's define a "chirp" signal:
 
%% Cell type:code id:72ac53c6 tags:
 
``` julia
function chirp(t, T, A, phi, omega, beta)
signal = A .* sin.(phi .+ omega .* (t .- T) .+ beta .* (t .- T).^2)
signal[t .< T] .= 0
return signal
end
 
N = 4096
t = range(-50,50,N)
h_true = chirp(t, -20, 0.8, 0, 0.2, 0.02)
d = Normal(0.,1.)
h = h_true .+ rand(d,N);
```
 
%% Cell type:code id:4b53eaa1 tags:
 
``` julia
fg = Figure()
 
ax1 = Axis(fg[1, 1],
xlabel="Time (s)",
ylabel="Flux (arbitrary units)",
)
 
lines!(t,h,label="Chirp signal")
lines!(t,h_true,label="Chirp signal without noise")
 
axislegend()
 
fg
```
 
%% Output
 
 
%% Cell type:markdown id:af4dc680 tags:
 
- Let's now compute a Fourier periodogram
 
%% Cell type:code id:c940d327 tags:
 
``` julia
# We use the DSP package
 
dt = 100/4096
psd = periodogram(h, fs=1/dt)
 
fg = Figure()
 
ax1 = Axis(fg[1, 1],
xlabel = "Frequency (Hz)",
ylabel = "Power",
)
 
lines!(psd.freq,psd.power,color=:blue,label="PSD")
 
xlims!(0,1)
 
axislegend()
 
 
fg
```
 
%% Output
 
 
%% Cell type:markdown id:f448373b tags:
 
- The periodogram shows power at the lowest frequencies with a break at frequencies higher than $\sim 0.5$Hz, i.e. for time-scales shorter than $\sim 1$s.
 
- Let's try to study the time-evolution of the PSD by a *spectrogram*:
 
%% Cell type:code id:b7ef3b6c tags:
 
``` julia
tper = spectrogram(h,div(length(h),5),0,fs=1/dt)
 
fg = Figure()
 
ax = Axis(fg[1, 1],
xlabel="Time (s)",
ylabel="Frequency (Hz)",
)
 
p = heatmap!(tper.time.-50,tper.freq,tper.power',colorscale=log10)
Colorbar(fg[:, end+1], p)
 
ylims!(0,1)
 
fg
```
 
%% Output
 
 
%% Cell type:markdown id:2455afd2 tags:
 
- And, it is indeed clearly visible a signal appearing at about -20s with frequncy increasing with time.
 
- Nevertheless, this is not an optimal solution since the number of bins is arbitrary, and clearly we reduce the resolution of the analysis decreasing the length of the time-series.
 
<br>
 
- This is where a wavelet analysis can give a contribution.
 
%% Cell type:code id:2b5fa16f tags:
 
``` julia
c = wavelet(Morlet(π), Q=8, β=1, averagingType=NoAve(), p=4)
 
res = cwt(h, c);
```
 
%% Cell type:code id:2186d3e3 tags:
 
``` julia
fg = Figure()
 
ax = Axis(fg[1, 1],
xlabel=L"Time index $\rightarrow$",
ylabel=L"Frequency index $\rightarrow$",
xticklabelsvisible=false,
yticklabelsvisible=false,
)
 
p = heatmap!(abs.(res))
Colorbar(fg[:, end+1], p)
 
fg
```
 
%% Output
 
 
%% Cell type:markdown id:87ff466d tags:
 
- And we can appreciate how the formation and evolution of the signal is clealry visible in the plot.
 
%% Cell type:markdown id:ffe137ba tags:
 
## Reference & Material
 
Material and papers related to the topics discussed in this lecture.
 
- [Belloni & Bhattacharya (2022) - "Basics of Fourier Analysis for High-Energy Astronomy”](https://ui.adsabs.harvard.edu/abs/2022hxga.book....7B/abstract)
- [Ivezić et al. (2020) - "Statistics, Data Mining, and Machine Learning in Astronomy"](https://ui.adsabs.harvard.edu/abs/2020sdmm.book.....I/abstract)
 
%% Cell type:markdown id:dff2c952 tags:
 
## Further Material
 
Papers for examining more closely some of the discussed topics.
 
- [Torrence & Compo (1998) - "A Practical Guide to Wavelet Analysis"](https://ui.adsabs.harvard.edu/abs/1998BAMS...79...61T/abstract)
 
 
%% Cell type:markdown id:2a98318f tags:
 
### Credits
***
 
This notebook contains material obtained from https://www.astroml.org/book_figures/chapter10/fig_chirp2_PSD.html.
 
%% Cell type:markdown id:05e93b1d tags:
 
## Course Flow
 
<table>
<tr>
<td>Previous lecture</td>
<td>Next lecture</td>
</tr>
<tr>
<td><a href="../Science%20Case%20-%20AGN%20and%20Blazars/Lecture-AGN-and-Blazars.ipynb">Science case about AGN and blazars</a></td>
<td><a href="../Lecture%20-%20Time%20of%20Arrival/Lecture-Time-of-Arrival.ipynb">Lecture about time of arrival analysis</a></td>
</tr>
</table>
 
 
%% Cell type:markdown id:591bd355 tags:
 
**Copyright**
 
This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the examples, unless obtained from other properly quoted sources, under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Stefano Covino, Time Domain Astrophysics - Lecture notes featuring computational examples, 2024*.
This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the examples, unless obtained from other properly quoted sources, under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Stefano Covino, Time Domain Astrophysics - Lecture notes featuring computational examples, 2025*.
%% Cell type:code id:14e835e0-2309-45de-b1cb-043de756eb48 tags:
``` julia
```
......
### A Pluto.jl notebook ###
# v0.20.4
using Markdown
using InteractiveUtils
# ╔═╡ 05ad5ab7-70f8-4bdb-8bbf-e05822cc5623
# ╠═╡ show_logs = false
import Pkg; Pkg.activate(".")
# ╔═╡ d2361176-db34-450c-8721-ed4391b90e9c
begin
using CairoMakie
using CommonMark
using ContinuousWavelets
using Distributions
using DSP
using PlutoUI
using Random
end
# ╔═╡ 0118270a-090b-4237-9002-bd68b608a307
md"""
**What is this?**
*This jupyter notebook is part of a collection of notebooks on various topics discussed during the Time Domain Astrophysics course delivered by Stefano Covino at the [Università dell'Insubria](https://www.uninsubria.eu/) in Como (Italy). Please direct questions and suggestions to [stefano.covino@inaf.it](mailto:stefano.covino@inaf.it).*
"""
# ╔═╡ de92ed35-d7b6-40f9-a642-d9a819804c85
md"""
**This is a `Julia` notebook**
"""
# ╔═╡ 1f00bd93-c681-4b4f-ba15-f4f67b3b40e3
Pkg.instantiate()
# ╔═╡ 5049cda4-f875-4291-a7e5-59a09b57e4e8
# ╠═╡ show_logs = false
Pkg.instantiate()
# ╔═╡ f40ad72a-b6fd-4024-ab09-de1b6e31501b
md"""
$(LocalResource("Pics/TimeDomainBanner.jpg"))
"""
# ╔═╡ 090ffe2f-b381-49ec-bc08-53219c690e77
md"""
# Wavelet Analysis
***
## PSD Varying Time-Series
***
- The trigonometric basis functions used in the Fourier transform have an infinite extent and for this reason the Fourier transform may not be the best method to analyze non periodic time series data, such as the case of a localized event (e.g., a burst that decays over some time scale so that the PSD is also varying with time).
- We can (and actually do) evaluate the PSD for finite stretches of time series in order to detect its changes.
- This approach (called *spectrogram*, or *dynamical power spectra analysis*, or even *windowed Fourier transform*) suffers from degraded spectral resolution and is sensitive to the specific choice of time series segmentation length.
### Windowed Fourier transform
***
- The FT is performed on a sliding segment of length $T$ from a time series of time step $δt$ and total length $Nδt$, thus returning frequencies from $T^{−1}$ to $2δt^{−1}$ at each time step.
- The segments can be windowed with an arbitrary function such as a boxcar (no smoothing) or any other choice.
- The WFT represents an inaccurate and inefficient method of time–frequency localization, as it imposes a scale or “response interval” T into the analysis.
- The inaccuracy arises from the aliasing of high- and low-frequency components that do not fall within the frequency range of the window. The inefficiency comes from the $T/(2δt)$ frequencies, which must be analyzed at each time step.
"""
# ╔═╡ 6721746e-71e6-4d3d-a32a-6c0385986351
md"""
## Wavelets
***
- A popular family of basis functions is called wavelets.
- By construction, wavelets are localized in both frequency and time domains.
- Individual wavelets are specified by a set of *wavelet filter coefficients*. Given a wavelet, a complete orthonormal set of basis functions can be constructed by scalings and translations.
- Different wavelet families trade the localization of a wavelet with its smoothness. Popular wavelets include “Mexican hat”, Haar and Daubechies wavelets.
- A wavelet is a “small wave,” where small derives from the fact that it is mostly limited to an interval of time.
- A (mother) wavelet $ψ$ has to satisfy two requirements: its integral must be zero and the integral of its square must be unity:
```math
\int_{-\infty}^\infty \Psi(t) dt = 0 \qquad \int_{-\infty}^\infty \Psi^2(t) dt = 1
```
- It is easy to see that these requirements imply that $ψ(t)$ is essentially non-zero only over a limited range of t and that it has to extend both above and below zero.
$(LocalResource("Pics/wavelets.png"))
- Three examples of wavelets. From left to right: Haar wavelet, mexican hat and morlet wavelet.
"""
# ╔═╡ 0bac9591-9018-47af-b567-ec756f8c6c42
cm"""
- A wavelet can be shifted by ``τ`` and dilated by a scale parameter ``σ``:
```math
\Psi_{\tau,\sigma}(t) = \frac{1}{\sqrt{\sigma}} \Psi\left(\frac{t-\tau}{\sigma}\right)
```
- Again ``ψ(t)`` is essentially non-zero only over a limited range of ``t`` and that it has to extend both above and below zero.
- The wavelet transform of a function ``f(t)`` is computed by correlating ``f(t)`` with the complex conjugate of ``ψ_{τ,σ}(t)`` (wavelets can also be complex functions):
```math
W[f(\tau,\sigma)] = \int_{-\infty}^\infty f(t) \Psi_{\tau,\sigma}^*(t)dt
```
- The *discrete wavelet transform* (DWT) can be used to analyze the power spectrum of a time series as a function of time.
- A possible wavelet might be: ``w(t|t_0, f_0, Q) = A \exp[i 2 \pi f_0 (t-t_0)] \exp[-f_0^2 (t-t_0)^2 / Q^2]``,
- where ``t_0`` is the central time, ``f_0`` is the central frequency, and the dimensionless parameter ``Q`` is a model parameter which controls the width of the frequency window.
$(LocalResource("Pics/waveletex1.png"))
- The Fourier transform of the previous wavelet is:
```math
W(t|t_0, f_0, Q) = \left( \frac{\pi}{f_0^2/Q^2} \right)^{1/2} \exp[i 2 \pi f_0 t_0] \exp[-\frac{\pi^2 Q^2 (f-f_0)^2}{Q f_0^2}]
```
- Technically speaking, this is not exactly a wavelet. It is indeed closely related to a true wavelet, the Morlet wavelet, through a simple scaling and offset.
- One might think to it as a “matched filters”.
- Given a signal ``h(t)``, its wavelet transform is given by the convolution:
```math
H_w(t_0; f_0, Q) = \int_{-\infty}^\infty h(t)w(t|t_0, f_0, Q)
```
- By the convolution theorem, we can write the Fourier transform of ``H_w`` as the pointwise product of the Fourier transforms of ``h(t)`` (e.g. by a DFT) and of ``w(t | t_0, f_0, Q)``.
- Then, the wavelet PSD, can be defined by:
```math
{\rm PSD}_w(f_0, t_0; Q) = |H_w(t_0; f_0, Q)|^2
```
- As said, unlike the typical Fourier-transform PSD, the wavelet PSD allows detection of frequency information which is localized in time.
- This is indeed one the approaches (we are simplifying a lot here…) used in the, e.g., LIGO/Virgo projects to detect gravitational wave events.
- Because of the noise level in the GW measurements, rather than a standard wavelet they instead use functions which are tuned to the expected form of the signal (i.e., matched filters).
- Example with a localized Gaussian noise:
$(LocalResource("Pics/waveletex2.png"))
- The upper panel shows the signal. The middle panel shows an example wavelet and the lower panel shows the power spectral density as a function of the frequency ``f_0`` and the time ``t_0``, for ``Q = 1.0``.
- And now a different example:
(LocalResource("Pics/waveletex3.png"))
- The upper panel shows the input signal, which consists of a Gaussian spike in the presence of white (Gaussian) noise. The middle panel shows again an example wavelet. The lower panel shows the PSD as a function of the frequency ``f_0`` and the time ``t_0``, for ``Q = 0.3``.
"""
# ╔═╡ a37c5419-797a-4626-8138-2f2d5ab8f416
md"""
#### Exercize: wavelet analysis of a chirp signal
***
- Let's define a "chirp" signal:
"""
# ╔═╡ f3f93bff-1c39-4652-aa41-0b7d88ad9caf
begin
function chirp(t, T, A, phi, omega, beta)
signal = A .* sin.(phi .+ omega .* (t .- T) .+ beta .* (t .- T).^2)
signal[t .< T] .= 0
return signal
end
N = 4096
t = range(-50,50,N)
h_true = chirp(t, -20, 0.8, 0, 0.2, 0.02)
d = Normal(0.,1.)
h = h_true .+ rand(d,N)
end;
# ╔═╡ 13fbe8ec-5266-4b72-b85e-f171c8c7ba90
begin
fg1 = Figure()
ax1fg1 = Axis(fg1[1, 1],
xlabel="Time (s)",
ylabel="Flux (arbitrary units)",
)
lines!(t,h,label="Chirp signal")
lines!(t,h_true,label="Chirp signal without noise")
axislegend()
fg1
end
# ╔═╡ 8550f405-2819-4a33-9cd2-c72dbbb04a87
md"""
- Let's now compute a Fourier periodogram
"""
# ╔═╡ 0e3b9aaf-0d91-4c60-8118-5624273d3983
begin
# We use the DSP package
dt = 100/4096
psd = periodogram(h, fs=1/dt)
fg2 = Figure()
ax1fg2 = Axis(fg2[1, 1],
xlabel = "Frequency (Hz)",
ylabel = "Power",
)
lines!(psd.freq,psd.power,color=:blue,label="PSD")
xlims!(0,1)
axislegend()
fg2
end
# ╔═╡ 7b3e4dfa-8dcd-4b5c-bc05-e1d51449d8c3
md"""
- The periodogram shows power at the lowest frequencies with a break at frequencies higher than $\sim 0.5$Hz, i.e. for time-scales shorter than $\sim 1$s.
- Let's try to study the time-evolution of the PSD by a *spectrogram*:
"""
# ╔═╡ 738b344d-e2ae-4f99-b755-38c0dcdc3da3
begin
tper = spectrogram(h,div(length(h),5),0,fs=1/dt)
fg3 = Figure()
axfg3 = Axis(fg3[1, 1],
xlabel="Time (s)",
ylabel="Frequency (Hz)",
)
p = heatmap!(tper.time.-50,tper.freq,tper.power',colorscale=log10)
Colorbar(fg3[:, end+1], p)
ylims!(0,1)
fg3
end
# ╔═╡ 576eb9f0-40c6-4a9a-9504-d562452541e6
md"""
- And, it is indeed clearly visible a signal appearing at about -20s with frequncy increasing with time.
- Nevertheless, this is not an optimal solution since the number of bins is arbitrary, and clearly we reduce the resolution of the analysis decreasing the length of the time-series.
- This is where a wavelet analysis can give a contribution.
"""
# ╔═╡ 6e81c233-8aa5-4ab4-a011-c56ad0026ea9
begin
c = wavelet(Morlet(π), Q=8, β=1, averagingType=NoAve(), p=4)
res = cwt(h, c)
end;
# ╔═╡ 6a71ccd0-8ee7-4c39-b8c7-88b398e98d1a
begin
fg4 = Figure()
axfg4 = Axis(fg4[1, 1],
xlabel=L"Time index $\rightarrow$",
ylabel=L"Frequency index $\rightarrow$",
xticklabelsvisible=false,
yticklabelsvisible=false,
)
p4 = heatmap!(abs.(res))
Colorbar(fg4[:, end+1], p4)
fg4
end
# ╔═╡ f153fdfa-57d5-4112-8d59-847ee1ec06ca
md"""
- And we can appreciate how the formation and evolution of the signal is clealry visible in the plot.
"""
# ╔═╡ bde447ae-3998-4a1a-8e4e-32f97b2854a1
md"""
## Reference & Material
Material and papers related to the topics discussed in this lecture.
- [Belloni & Bhattacharya (2022) - "Basics of Fourier Analysis for High-Energy Astronomy”](https://ui.adsabs.harvard.edu/abs/2022hxga.book....7B/abstract)
- [Ivezić et al. (2020) - "Statistics, Data Mining, and Machine Learning in Astronomy"](https://ui.adsabs.harvard.edu/abs/2020sdmm.book.....I/abstract)
"""
# ╔═╡ c4b73630-ba8d-4403-9b7b-b498eac9792b
md"""
## Further Material
Papers for examining more closely some of the discussed topics.
- [Torrence & Compo (1998) - "A Practical Guide to Wavelet Analysis"](https://ui.adsabs.harvard.edu/abs/1998BAMS...79...61T/abstract)
"""
# ╔═╡ af38d725-6e0f-4958-85bb-1be029ba3fb5
md"""
### Credits
***
This notebook contains material obtained from [https://www.astroml.org/book_figures/chapter10/fig_chirp2_PSD.html](https://www.astroml.org/book_figures/chapter10/fig_chirp2_PSD.html).
"""
# ╔═╡ 241907cb-1fb0-4963-8bb9-7894136e43de
cm"""
## Course Flow
<table>
<tr>
<td>Previous lecture</td>
<td>Next lecture</td>
</tr>
<tr>
<td><a href="./open?path=Lectures/Science Case - AGN and Blazars/Lecture-AGN-and-Blazars.jl">Science case about AGN and blazars</a></td>
<td><a href="./open?path=Lecture - Time of Arrival/Lecture-Time-of-Arrival.jl">Lecture about time of arrival analysis</a></td>
</tr>
</table>
"""
# ╔═╡ d3b81af7-b85b-4212-82a2-0f97ec31a484
md"""
**Copyright**
This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the examples, unless obtained from other properly quoted sources, under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Stefano Covino, Time Domain Astrophysics - Lecture notes featuring computational examples, 2025*.
"""
# ╔═╡ Cell order:
# ╟─0118270a-090b-4237-9002-bd68b608a307
# ╟─de92ed35-d7b6-40f9-a642-d9a819804c85
# ╟─05ad5ab7-70f8-4bdb-8bbf-e05822cc5623
# ╟─1f00bd93-c681-4b4f-ba15-f4f67b3b40e3
# ╟─5049cda4-f875-4291-a7e5-59a09b57e4e8
# ╟─d2361176-db34-450c-8721-ed4391b90e9c
# ╟─f40ad72a-b6fd-4024-ab09-de1b6e31501b
# ╟─090ffe2f-b381-49ec-bc08-53219c690e77
# ╟─6721746e-71e6-4d3d-a32a-6c0385986351
# ╟─0bac9591-9018-47af-b567-ec756f8c6c42
# ╟─a37c5419-797a-4626-8138-2f2d5ab8f416
# ╟─f3f93bff-1c39-4652-aa41-0b7d88ad9caf
# ╟─13fbe8ec-5266-4b72-b85e-f171c8c7ba90
# ╟─8550f405-2819-4a33-9cd2-c72dbbb04a87
# ╟─0e3b9aaf-0d91-4c60-8118-5624273d3983
# ╟─7b3e4dfa-8dcd-4b5c-bc05-e1d51449d8c3
# ╟─738b344d-e2ae-4f99-b755-38c0dcdc3da3
# ╟─576eb9f0-40c6-4a9a-9504-d562452541e6
# ╟─6e81c233-8aa5-4ab4-a011-c56ad0026ea9
# ╟─6a71ccd0-8ee7-4c39-b8c7-88b398e98d1a
# ╟─f153fdfa-57d5-4112-8d59-847ee1ec06ca
# ╟─bde447ae-3998-4a1a-8e4e-32f97b2854a1
# ╟─c4b73630-ba8d-4403-9b7b-b498eac9792b
# ╟─af38d725-6e0f-4958-85bb-1be029ba3fb5
# ╟─241907cb-1fb0-4963-8bb9-7894136e43de
# ╟─d3b81af7-b85b-4212-82a2-0f97ec31a484
This diff is collapsed.
[deps]
CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
CommonMark = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"
ContinuousWavelets = "96eb917e-2868-4417-9cb6-27e7ff17528f"
DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
......@@ -2375,7 +2375,7 @@
"source": [
"**Copyright**\n",
"\n",
"This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the examples, unless obtained from other properly quoted sources, under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Stefano Covino, Time Domain Astrophysics - Lecture notes featuring computational examples, 2024*."
"This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources). Feel free to use the notebook for your own purposes. The text is licensed under [Creative Commons Attribution 4.0](https://creativecommons.org/licenses/by/4.0/), the code of the examples, unless obtained from other properly quoted sources, under the [MIT license](https://opensource.org/licenses/MIT). Please attribute the work as follows: *Stefano Covino, Time Domain Astrophysics - Lecture notes featuring computational examples, 2025*."
]
},
{
......@@ -2392,7 +2392,7 @@
"provenance": []
},
"kernelspec": {
"display_name": "Julia 1.11.3",
"display_name": "Julia 1.11.4",
"language": "julia",
"name": "julia-1.11"
},
......@@ -2400,7 +2400,7 @@
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.11.3"
"version": "1.11.4"
}
},
"nbformat": 4,
......@@ -2,4 +2,4 @@
This is a repository with material (notebooks, papers, etc.) for the **Time Domain Astrophysics** course delivered at the *Università dell'Insubria* by Stefano Covino.
*Last update: 11 March 2025.*
*Last update: 12 March 2025.*
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment