Created
January 4, 2021 23:36
-
-
Save vly/73f7910a68ac6ca3ec1ee3fc0eccaca3 to your computer and use it in GitHub Desktop.
Henderson filter in Julia
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
""" | |
Implementation of Henderson filter in Julia. | |
Followed https://github.com/bpalmer4 python implementation | |
""" | |
using DataFrames | |
function hmaSymmetricWeights(n::Int) | |
m = (n-1)÷2 | |
m1 = (m+1)^2 | |
m2 = (m+2)^2 | |
m3 = (m+3)^2 | |
d = float(8*(m+2)*(m2-1)*(4*m2-1)*(4*m2-9)*(4*m2-25)) | |
w = map(x -> (315*(m1-(x^2))*(m2-(x^2))*(m3-(x^2))*(3*m2-11*(x^2)-16))/d, range(0, m+1, step = 1))[end-1:-1:1] | |
u = vcat(w, w[end-1:-1:1]) | |
return if n % 2 > 0 | |
u | |
else | |
vcat(u, missing) | |
end | |
end | |
function hmaAsymmetricWeights(m::Int, w::AbstractArray{<:Number}) | |
n = length(w) | |
@assert m <= n "The m argument must be less than n" | |
@assert m >= (n-1)÷2 "The m argument must be greater than (n-1)/2" | |
sumResidual = sum(w[range(m + 1, n, step = 1)]) / m | |
sumEnd = sum(map(x -> (x - ((m+1.0)/2.0)) * w[x], range(m+1, n, step = 1))) | |
ic = if n >= 13 && n < 15 | |
3.5 | |
elseif n >= 15 | |
4.5 | |
else 1.0 | |
end | |
b2s2 = 4/pi/ic^2 | |
denominator = 1.0 + ((m*(m-1.0)*(m+1.0) / 12.0 ) * b2s2) | |
u = map(r -> w[r] + sumResidual + (((r-1+1.0) - (m+1.0)/2.0) * b2s2) / denominator * sumEnd, | |
range(1, m, step = 1)) | |
return (u) | |
end | |
""" | |
Henderson(s, n) | |
Applies the Henderson filter to dataset `s` with `n` term. | |
""" | |
function Henderson(s::AbstractArray{<:Number}, n::Int) | |
@assert isodd(n) "n must be odd" | |
@assert n >= 5 "n must be >= 5" | |
@assert length(s) >= n "dataset must be >= than n" | |
w = hmaSymmetricWeights(n) | |
m = (n-1) ÷ 2 | |
l = length(s) | |
map(i -> if i - 1 < m | |
u = hmaAsymmetricWeights(m + i, w)[end:-1:1] | |
sum(s[1:i + m] .* u) | |
elseif i - 1 + m >= l | |
u = hmaAsymmetricWeights(m + l - i + 1, w) | |
sum(s[i-m:l] .* u) | |
else | |
sum(s[i-m:i+m] .* w) | |
end, range(1, l, step = 1)) | |
end |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment