diff --git a/LICENSE b/LICENSE index 4c469e2..f60126c 100644 --- a/LICENSE +++ b/LICENSE @@ -1,6 +1,6 @@ MIT License -Copyright (c) 2022 Ujjwal Panda +Copyright (c) 2022-2023 Ujjwal Panda 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 index 9e9ba27..3d94362 100644 --- 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 "] @@ -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 index d0a45f3..8b13789 100644 --- a/README.md +++ b/README.md @@ -1,39 +1 @@ -
-

- -Noisy : -Simulating power law noise in Julia. - -

-
- -![License][license] -![GitHub Stars][stars] -[![Gitmoji Badge][gitmoji_badge]][gitmoji] - -
- -
- -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) - -
-
- -[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 deleted file mode 100644 index 2dd3376..0000000 Binary files a/noisy.png and /dev/null differ diff --git a/noisy2d.png b/noisy2d.png deleted file mode 100644 index 062bd35..0000000 Binary files a/noisy2d.png and /dev/null differ diff --git a/pix/noisy1D.png b/pix/noisy1D.png new file mode 100644 index 0000000..92b1d5b Binary files /dev/null and b/pix/noisy1D.png differ diff --git a/pix/noisy2D.png b/pix/noisy2D.png new file mode 100644 index 0000000..460f3e0 Binary files /dev/null and b/pix/noisy2D.png differ diff --git a/pix/spectrum.png b/pix/spectrum.png new file mode 100644 index 0000000..d51d88f Binary files /dev/null and b/pix/spectrum.png differ diff --git a/src/Noisy.jl b/src/Noisy.jl deleted file mode 100644 index c575097..0000000 --- a/src/Noisy.jl +++ /dev/null @@ -1,51 +0,0 @@ -module Noisy - -export noise - -using FFTW -using Random - -zerofirst!(A) = selectdim(A, ndims(A), lastindex(A, ndims(A))) .= 0.0 -zerolast!(A) = selectdim(A, ndims(A), firstindex(A, ndims(A))) .= 0.0 - -function scale!(A, f, f₀, β) - for I ∈ LinearIndices(A) - A[I] = (f[I] < f₀ ? f₀ : f[I])^(-β/2) - end -end - -function σ(A, nf, nt) - Σ = zero(eltype(A)) - for I ∈ LinearIndices(A) - Σ += 1 < I < nf ? A[I]^2 : 0.0 - end - Σ += (A[end] * (1 + (nt % 2)) / 2)^2 - 2 * √Σ / nt -end - -function noise!(X, A, nt) - x = randn(Float64, size(X)...) - y = randn(Float64, size(X)...) - (nt % 2) == 0 && zerolast!(y) - zerofirst!(y) - for I in CartesianIndices(X) - X[I] = A[Tuple(I)[end]] * (x[I] + (im * y[I])) - end - irfft(X, nt, ndims(X)) ./ σ(A, size(A, 1), nt) -end - -function noise(β, f₀, dims...) - nt = dims[end] - f = rfftfreq(nt) - f₀ = 0 <= f₀ <= 0.5 ? max(f₀, 1/nt) : throw("ERROR: f₀ ∉ [0, 0.5].") - X = Array{ComplexF64}(undef, dims[1:end-1]..., size(f, 1)) - A = similar(f) - scale!(A, f, f₀, β) - noise!(X, A, nt) -end - -""" -""" -noise - -end diff --git a/src/PowerLawNoise.jl b/src/PowerLawNoise.jl new file mode 100644 index 0000000..9ee847b --- /dev/null +++ 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