Commit 6fadd8cb authored by Stefano Covino's avatar Stefano Covino
Browse files

110325 commit

parent 8b834d90
Loading
Loading
Loading
Loading

Course backup 1.jl

0 → 100644
+0 −0

Empty file added.

+4 −4
Original line number Diff line number Diff line
@@ -57,9 +57,9 @@ md"""
5. [Science case: X-ray binaries](./open?path=Lectures/Science Case - X-Ray Binaries/Lecture-X-RayBinaries.jl)
6. [Lecture: Irregular sampling](./open?path=Lectures/Lecture - Lomb-Scargle/Lecture-Lomb-Scargle.jl)
7. [Science case: Variable stars](./open?path=Lectures/Science Case - Variable Stars/Lecture-VariableStars.jl)
8. [Lecture: Time Domain analysis](./open?path=Lectures/Lecture - Time Domain%20Analysis/Lecture-Time-Domain.jl)
9. [Science case: AGN and Blazars](Lectures/Science%20Case%20-%20AGN%20and%20Blazars/Lecture-AGN-and-Blazars.ipynb)
10. [Lecture: Wavelet Analysis](Lectures/Lecture%20-%20Wavelet%20Analysis/Lecture-Wavelet-Analysis.ipynb)
8. [Lecture: Time Domain analysis](./open?path=Lectures/Lecture - Time Domain Analysis/Lecture-Time-Domain.jl)
9. [Science case: AGN and Blazars](./open?path=Lectures/Science Case - AGN and Blazars/Lecture-AGN-and-Blazars.jl)
10. [Lecture: Wavelet Analysis](./open?path=Lectures/Lecture - Wavelet Analysis/Lecture-Wavelet-Analysis.jl)
11. [Lecture: Time of Arrival](Lectures/Lecture%20-%20Time%20of%20Arrival/Lecture-Time-of-Arrival.ipynb)
12. [Science case: FRBs](Lectures/Science%20Case%20-%20FRBs/Lecture-FRBs.ipynb)
13. [Lecture: Non Parametric Analysis](Lectures/Lecture%20-%20Non%20Parametric%20Analysis/Lecture-NonParametricAnalysis.ipynb)
@@ -88,7 +88,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"

+21 −3
Original line number Diff line number Diff line
@@ -4,6 +4,18 @@
using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    #! format: off
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
    #! format: on
end

# ╔═╡ dfdaa2c2-58bc-4911-8e3b-de854e5d56fa
# ╠═╡ show_logs = false
import Pkg; Pkg.activate(".")
@@ -63,6 +75,11 @@ md"""
 - where we have defined the "sinc" function defined as ${\rm sinc} (t) = \sin(\pi t) / \pi t$.
"""

# ╔═╡ c7364a7d-6d80-4872-ab16-b07cbc659993
md"""
Width of the box function: $( @bind wdth PlutoUI.Slider(0.1:0.1:4, default=1) ) 
"""

# ╔═╡ fc298404-78c0-49b9-8250-091dbc2b7487
begin
	T=2.
@@ -86,8 +103,8 @@ begin
	    xlabel = "t",
	    )
	
	T=1
	lines!(xrng,map(x -> f(x;T=T),xrng),label="T = "*string(T))
	T=wdth
	lines!(xrng,map(x -> f(x;T=wdth),xrng),label="T = "*string(T))
	T=3
	lines!(xrng,map(x -> f(x;T=3),xrng),label="T = "*string(T))
	
@@ -102,7 +119,7 @@ begin
	    xlabel = "f",
	    )
	
	T=1
	T=wdth
	lines!(frng,T*sinc.(T*frng),label="T = "*string(T))
	T=3
	lines!(frng,T*sinc.(T*frng),label="T = "*string(T))
@@ -204,6 +221,7 @@ This notebook is provided as [Open Educational Resource](https://en.wikipedia.or
# ╟─d34ab604-27fa-4877-a8e1-d17bb8fb026f
# ╟─48232ae4-81b0-4cf4-80a9-1a0085ba8d9f
# ╟─e5b7c371-b30d-44d3-8ee9-f347719ba366
# ╟─c7364a7d-6d80-4872-ab16-b07cbc659993
# ╠═fc298404-78c0-49b9-8250-091dbc2b7487
# ╟─9b5f3662-37da-4e19-a4a3-a80dbbfb041a
# ╟─4a49d01f-0138-4ecd-b5af-a9328cfcbd5c
+77 −22
Original line number Diff line number Diff line
@@ -4,6 +4,18 @@
using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    #! format: off
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
    #! format: on
end

# ╔═╡ 419e9d2c-f5df-49a3-b0ce-9804aacda4b7
# ╠═╡ show_logs = false
import Pkg; Pkg.activate(".")
@@ -21,6 +33,7 @@ begin
	using Format
	using HTTP
	using LaTeXStrings
	using Latexify
	using PlutoUI
	using Statistics
end
@@ -53,12 +66,12 @@ cm"""
***

- A time-series is any sequene of observation such that the distribution of a given value depends on the previous values.
- Time is an exogeneous (outside the model) variable that is directional - measuremets only depend on the past.
- Time is an exogeneous (outside the model) variable that is directional - measurements only depend on the past.
    - This is a statement of causality.

- However, the exogenous variable can be anything.

- Let'assume to have a set of data extracted from ``y(t) = A \sin(\omega t)`` with homoscedastic variance ``V = \sigma^2 + A^2/2``.
- Let's assume to have a set of data extracted from ``y(t) = A \sin(\omega t)`` with homoscedastic variance ``V = \sigma^2 + A^2/2``.
    - This is easy to prove if you compute the variance as ``\sum (y-\lt y \gt)^2 / N``. Since the average value is zero, this turns out to be ``V = \frac{A^2}{N} \sum \sin^2 (\omega t)`` giving the ``A^2/2`` term.

- We can compute the ``\chi^2`` for this toy model:
@@ -138,7 +151,7 @@ f(t) = {1\over 2\pi} \int\limits_{-\infty}^{+\infty} F(\omega) e^{i\omega t} d\o

| $f(t)$ | F(ω) |
|:-------: |:--------:|
| Real    | H(-f) = H$^*$(ω)   |
| Real    | H(-ω) = H$^*$(ω)   |
| Even    | Even   |
| Odd    | Odd   |
| Real and Even | Real and Even     |
@@ -183,7 +196,11 @@ P(\omega) = F(\omega)\cdot F^*(\omega) = | F(\omega) |^2
- E.g., if two signals are $f(t)$ and $g(t)$, the PSD of the sum of the two is:

```math
P[f(t)+g(t)] = |F[f(t)+g(t)]|^2 =  F[f(t)+g(t)]\cdot F^*[f(t)+g(t)] = P[f(t)] + P[g(t)] + 2 Re\{F[f(t)]\cdot F[g(t)] \}
P[f(t)+g(t)] = |F[f(t)+g(t)]|^2 =  F[f(t)+g(t)]\cdot F^*[f(t)+g(t)] = 
```

```math
= P[f(t)] + P[g(t)] + 2 Re\{F[f(t)]\cdot F[g(t)] \}
```

- If the two signals are uncorrelated, the cross-term is zero and linearity applies.
@@ -351,9 +368,15 @@ md"""
begin
	posfr = temp_freq .> 0
	
	println("Nyquist frequency: ", round(maximum(temp_freq) / 365, digits=1), L" day$^{-1}$")
	#println("Nyquist frequency: ", round(maximum(temp_freq) / 365, digits=1), L" day$^{-1}$")
	
end
end;

# ╔═╡ 6f29671a-048a-49e1-b7f9-e5d9cd733a0c
Markdown.parse("""
##### Nyquist frequency:  $(latexify(maximum(temp_freq) / 365,fmt="%.1f"))

""")

# ╔═╡ b8a2074a-490e-4e92-ba59-3735184d7c91
md"""
@@ -381,6 +404,11 @@ md"""
- Or with a better zoom in a region of our interest.
"""

# ╔═╡ 0f5adef3-5de0-4fa3-b4ca-8ee59968047e
md"""
Plot horizontal axis limit: $( @bind xmx PlutoUI.Slider(0.1:0.5:200, default=7) ) 
"""

# ╔═╡ 1b13be55-5153-4d53-ba09-92832f67a448
begin
	fg3 = Figure()
@@ -392,7 +420,7 @@ begin
	
	lines!(temp_freq[posfr],10 * log10.(temp_psd[posfr]),color=:blue)
	
	xlims!(0,7)
	xlims!(0,xmx)
	ylims!(30,85)
	
	fg3
@@ -448,7 +476,7 @@ F[x \cdot y] = F[x] \otimes F[y] = \int\limits_{-\infty}^{+\infty}F[x(\nu')]F[y(
x_k = h(t) \cdot w(t) \cdot s(t)
```

- `w(t)` is a boxcar window function, which is 1 in the `(0,T)` interval and zero outside. `s(t)' is a series of delta functions at `t_k`, spaced by `T/N`:
- ``w(t)`` is a boxcar window function, which is 1 in the ``(0,T) `` interval and zero outside. ``s(t)`` is a series of delta functions at ``t_k``, spaced by ``T/N``:

$(LocalResource("Pics/windowing.png"))
"""
@@ -715,9 +743,18 @@ md"""
### Exercise about PSD manipulation
***

- Let's generate two arrays of relative timestamps, one 8 seconds long and one 1600 seconds long, with dt = 0.03125 s, and make two signals in units of counts. The signal is a sine wave with amplitude = 300 cts/s, frequency = 2 Hz, phase offset = 0 radians, and mean = 1000 cts/s. We then add Poisson noise to the light curve amd plot the shortest of the two.
- Let's generate two arrays of relative timestamps, one 8 seconds long and one 1600 seconds long, with dt = 0.03125 s, and make two signals in units of counts. The signal is a sine wave with amplitude = 300 cts/s, frequency = 2 Hz, phase offset = 0 radians, and mean = 1000 cts/s. We then add Poisson noise to the light curve and plot the shortest of the two.
"""

# ╔═╡ 1112b278-802d-4611-8954-2795db178e88
pf = @bind pl_fr NumberField(1:10, default=2);

# ╔═╡ 19a6c251-adc1-422d-a702-20ea9b971449
cm"Frequency for the test plot:"

# ╔═╡ a8b99222-18a7-4b47-8f98-8961a708ec57
pf

# ╔═╡ b96525bf-90e5-47b3-b305-1db0ad675c2c
begin
	dt = 0.03125  # seconds
@@ -726,8 +763,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.(2 .* pi .* times ./ 0.5) .+ 1000  # counts/s
	long_signal = 300 .* sin.(2 .* pi .* long_times ./ 0.5) .+ 1000  # counts/s
	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
	
	noisy = [rand(Poisson(theta)) for theta in signal .* dt]  # counts
	long_noisy = [rand(Poisson(theta)) for theta in long_signal .* dt]  # counts
@@ -752,6 +789,9 @@ md"""
- Now let's compute the periodogram
"""

# ╔═╡ b7213b38-21ac-4fb5-ab53-694d3ac2fda9
pf

# ╔═╡ 25fcf199-7da4-4632-a2fa-d0131f58d533
begin
	# We use the DSP package
@@ -785,11 +825,12 @@ md"""
- The zero frequency is the sum of the signal value, and it not typically used.
"""

# ╔═╡ 0770fef2-cade-46d5-9460-765e4fb0a1c9
begin
	println("Number of data points: ", length(noisy))
	println("Number of frequencies: ", length(psd.freq[2:end]))
end
# ╔═╡ 817993c0-a1ce-46d2-a213-ed6448d5159f
Markdown.parse("""
Number of data points:   $(latexify(length(noisy)))

Number of data points:   $(latexify(length(psd.freq[2:end])))
""")

# ╔═╡ 7fedf685-66f5-4d46-b52e-99c7305d3f95
md"""
@@ -798,6 +839,9 @@ md"""
- We want to average periodograms computed each 8 seconds, averaging therefore 1600/8 = 200 periodograms of 256 elements each.
"""

# ╔═╡ 9985012a-ad66-434b-84c2-13eead6a1af7
pf

# ╔═╡ 39538b6b-414d-436a-a1c4-19c6f6754dc7
begin
	apsd = welch_pgram(long_noisy, 256, 0, fs=1/dt)
@@ -826,6 +870,9 @@ md"""
- Let's try now to compute periodograms using different windows functions.
"""

# ╔═╡ cef2b9a7-701e-4975-8e71-f8a0afff1641
pf

# ╔═╡ b8c5baed-5a58-4380-b941-a7d4035c8d01
begin
	psd_rect = periodogram(noisy; fs=1/dt, window=nothing)
@@ -1157,9 +1204,11 @@ This notebook is provided as [Open Educational Resource](https://en.wikipedia.or
# ╠═b49f1cda-582a-4d43-9e0a-a0e6adf85e78
# ╟─c601e37c-c08d-4ce4-8843-6fd52bd1e812
# ╠═aff12d14-4803-480e-aab0-6d33917304b3
# ╟─6f29671a-048a-49e1-b7f9-e5d9cd733a0c
# ╟─b8a2074a-490e-4e92-ba59-3735184d7c91
# ╠═40fe7f7d-8d30-40f1-b6fa-932374c93380
# ╟─1eb98eee-8373-4fce-b40f-3118d096f082
# ╟─0f5adef3-5de0-4fa3-b4ca-8ee59968047e
# ╠═1b13be55-5153-4d53-ba09-92832f67a448
# ╟─bf2f9df2-7df7-43d2-926d-7a561be79c0e
# ╠═bea9ec97-57a7-45a4-9a29-ccca274fd5dc
@@ -1182,15 +1231,21 @@ This notebook is provided as [Open Educational Resource](https://en.wikipedia.or
# ╟─0b0b8578-4c67-463a-bd53-b3e6bd1f9712
# ╟─8e458457-f940-460a-8638-debe1f94defc
# ╟─b2e2bace-7e50-41f0-87fb-780737ff6616
# ╠═b96525bf-90e5-47b3-b305-1db0ad675c2c
# ╟─1112b278-802d-4611-8954-2795db178e88
# ╠═19a6c251-adc1-422d-a702-20ea9b971449
# ╟─a8b99222-18a7-4b47-8f98-8961a708ec57
# ╟─b96525bf-90e5-47b3-b305-1db0ad675c2c
# ╟─ae8a4ebf-da8d-4778-ad9a-87e3aab72399
# ╠═25fcf199-7da4-4632-a2fa-d0131f58d533
# ╟─b7213b38-21ac-4fb5-ab53-694d3ac2fda9
# ╟─25fcf199-7da4-4632-a2fa-d0131f58d533
# ╟─3e27b7d8-9dec-4e54-9050-5ca9275d8499
# ╠═0770fef2-cade-46d5-9460-765e4fb0a1c9
# ╟─817993c0-a1ce-46d2-a213-ed6448d5159f
# ╟─7fedf685-66f5-4d46-b52e-99c7305d3f95
# ╠═39538b6b-414d-436a-a1c4-19c6f6754dc7
# ╟─9985012a-ad66-434b-84c2-13eead6a1af7
# ╟─39538b6b-414d-436a-a1c4-19c6f6754dc7
# ╟─549cc9b7-01ff-4639-b8cb-66df220279a9
# ╠═b8c5baed-5a58-4380-b941-a7d4035c8d01
# ╟─cef2b9a7-701e-4975-8e71-f8a0afff1641
# ╟─b8c5baed-5a58-4380-b941-a7d4035c8d01
# ╟─ee646903-dc54-46d0-b2f5-b30e41a68853
# ╟─454fc15d-ea34-4938-b86e-76f2cb0d5125
# ╟─6e178c7b-1846-4390-b842-b4d303688bd8
+17 −1
Original line number Diff line number Diff line
@@ -2,7 +2,7 @@

julia_version = "1.11.3"
manifest_format = "2.0"
project_hash = "ec9b03d4059354e9185c514cbcd85f704a9629cc"
project_hash = "f29ba2d9e3fb794cb0b28e15f2cd85efe4c8b504"

[[deps.AbstractFFTs]]
deps = ["LinearAlgebra"]
@@ -750,6 +750,22 @@ git-tree-sha1 = "dda21b8cbd6a6c40d9d02a73230f9d70fed6918c"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.4.0"

[[deps.Latexify]]
deps = ["Format", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "OrderedCollections", "Requires"]
git-tree-sha1 = "cd714447457c660382fe634710fb56eb255ee42e"
uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
version = "0.16.6"

    [deps.Latexify.extensions]
    DataFramesExt = "DataFrames"
    SparseArraysExt = "SparseArrays"
    SymEngineExt = "SymEngine"

    [deps.Latexify.weakdeps]
    DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
    SymEngine = "123dc426-2d89-5057-bbad-38513e3affd8"

[[deps.LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
Loading