Skip to content

Instantly share code, notes, and snippets.

@vly
Created January 4, 2021 23:36
Show Gist options
  • Save vly/73f7910a68ac6ca3ec1ee3fc0eccaca3 to your computer and use it in GitHub Desktop.
Save vly/73f7910a68ac6ca3ec1ee3fc0eccaca3 to your computer and use it in GitHub Desktop.
Henderson filter in Julia
"""
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