🎨 Improve code structure etc.

astrogewgaw · May 2, 2023 · e4e15c0 · e4e15c0
1 parent b499b9b
commit e4e15c0
Show file tree

Hide file tree

Showing 10 changed files with 51 additions and 92 deletions.
diff --git a/LICENSE b/LICENSE
@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2022 Ujjwal Panda <[email protected]>
+Copyright (c) 2022-2023 Ujjwal Panda <[email protected]>
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

diff --git a/Project.toml b/Project.toml
@@ -1,4 +1,4 @@
-name = "Noisy"
+name = "PowerLawNoise"
 version = "0.0.1"
 uuid = "9abb9eec-80c7-4a18-accd-d591357fead6"
 authors = ["Ujjwal Panda <[email protected]>"]
@@ -8,4 +8,4 @@ FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
-julia = "1.7"
+julia = "1.8"
diff --git a/README.md b/README.md
@@ -1,39 +1 @@
-<div align="left" style="font-family:juliamono">
-<h2>
-<i>
-<code>Noisy</code> :
-<small><u>Simulating power law noise in Julia.</u></small>
-</i>
-</h2>
 
-<br/>
-
-![License][license]
-![GitHub Stars][stars]
-[![Gitmoji Badge][gitmoji_badge]][gitmoji]
-
-<br/>
-
-<div align="justify">
-
-This package implements a novel algorithm, following the work in [**"*On
-generating power law noise*", Timmer, J. and Koenig, M., Astron. Astrophys. 300,
-707-710 (1995)**][paper], which allows generating power law noise with arbitrary
-dimensions and arbitrary exponents. An exponent of 0, 1 and 2 corresponds to
-white, pink and red noise. The figure below demonstrates the difference between
-them, by plotting the resultant time series, generated via `Noisy.jl`:
-
-![Plot: 1D Noise](./noisy.png)
-
-Here is another one, but this time we are simulating noise in two dimensions:
-
-![Plot: 2D Noise](./noisy2d.png)
-
-</div>
-</div>
-
-[gitmoji]: https://gitmoji.dev
-[paper]: https://ui.adsabis.harvard.edu/abs/1995A%26A...300..707T/abstract
-[stars]: https://img.shields.io/github/stars/astrogewgaw/Noisy.jl?style=for-the-badge
-[license]: https://img.shields.io/github/license/astrogewgaw/Noisy.jl?style=for-the-badge
-[gitmoji_badge]: https://img.shields.io/badge/gitmoji-%20😜%20😍-FFDD67.svg?style=for-the-badge
diff --git a/noisy.png b/noisy.png
diff --git a/noisy2d.png b/noisy2d.png
diff --git a/pix/noisy1D.png b/pix/noisy1D.png
diff --git a/pix/noisy2D.png b/pix/noisy2D.png
diff --git a/pix/spectrum.png b/pix/spectrum.png
diff --git a/src/Noisy.jl b/src/Noisy.jl
diff --git a/src/PowerLawNoise.jl b/src/PowerLawNoise.jl
@@ -0,0 +1,48 @@
+module Noisy
+
+export noisy
+
+using FFTW
+using Random
+
+zerolast!(A) = selectdim(A, ndims(A), lastindex(A, ndims(A))) .= 0.0
+zerofirst!(A) = selectdim(A, ndims(A), firstindex(A, ndims(A))) .= 0.0
+
+function σ(A, n₀, n₁)
+    2 * √(
+        sum([1 < i < n₁ ? A[i]^2 : 0.0 for i ∈ eachindex(A)]) + (A[end] * (1 + (n₀ % 2)) / 2)^2,
+    ) / n₀
+end
+
+function scales!(A, ν, ν₀, β)
+    for i ∈ LinearIndices(A)
+        A[i] = (ν[i] < ν₀ ? ν₀ : ν[i])^(-β / 2)
+    end
+    A
+end
+
+function scales(ν, ν₀, β)
+    A = similar(ν)
+    scales!(A, ν, ν₀, β)
+end
+
+function noisy(β, ν₀, dims...)
+    n₀ = dims[end]
+    ν = rfftfreq(n₀)
+    ν₀ = 0 <= ν₀ <= 0.5 ? max(ν₀, n₀^-1) : throw("ν₀ ∉ [0, 0.5].")
+
+    N = Array{ComplexF64}(undef, dims[1:end-1]..., size(ν, 1))
+
+    X = randn(Float64, size(N)...)
+    Y = randn(Float64, size(N)...)
+
+    zerofirst!(Y)
+    A = scales(ν, ν₀, β)
+    (n₀ % 2) == 0 && zerolast!(Y)
+    for i in CartesianIndices(N)
+        N[i] = A[Tuple(i)[end]] * (X[i] + (im * Y[i]))
+    end
+    irfft(N, n₀, ndims(N)) ./ σ(A, size(A, 1), n₀)
+end
+
+end