Last active
April 26, 2024 12:03
-
-
Save jiahao/b8b5ac328c18b7ae8a17 to your computer and use it in GitHub Desktop.
An implementation of the Savitzky-Golay filter using generated functions. Accompanies https://medium.com/@acidflask/smoothing-data-with-julia-s-generated-functions-c80e240e05f3
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
""" | |
Savitzky-Golay filter of window half-width M and degree N | |
M is the number of points before and after to interpolate, i.e. the full width | |
of the window is 2M+1 | |
""" | |
immutable SavitzkyGolayFilter{M,N} end | |
@generated function Base.call{M,N,T}(::Type{SavitzkyGolayFilter{M,N}}, | |
data::AbstractVector{T}) | |
#Create Jacobian matrix | |
J = zeros(2M+1, N+1) | |
for i=1:2M+1, j=1:N+1 | |
J[i, j] = (i-M-1)^(j-1) | |
end | |
e₁ = zeros(N+1) | |
e₁[1] = 1.0 | |
#Compute filter coefficients | |
C = J' \ e₁ | |
#Evaluate filter on data matrix | |
To = typeof(C[1] * one(T)) #Calculate type of output | |
expr = quote | |
n = size(data, 1) | |
smoothed = zeros($To, n) | |
@inbounds for i in eachindex(smoothed) | |
smoothed[i] += $(C[M+1])*data[i] | |
end | |
smoothed | |
end | |
for j=1:M | |
insert!(expr.args[6].args[2].args[2].args, 1, | |
:(if i - $j ≥ 1 | |
smoothed[i] += $(C[M+1-j])*data[i-$j] | |
end) | |
) | |
push!(expr.args[6].args[2].args[2].args, | |
:(if i + $j ≤ n | |
smoothed[i] += $(C[M+1+j])*data[i+$j] | |
end) | |
) | |
end | |
return expr | |
end | |
#Two sample invocations | |
@show SavitzkyGolayFilter{2,3}([1:11;]) | |
@show SavitzkyGolayFilter{1,1}([1:11;]) |
I adapted it to a working update for Julia v1.1:
struct SavitzkyGolayFilter{M,N} end
@generated function (::SavitzkyGolayFilter{M,N})(data::AbstractVector{T}) where {M, N, T}
#Create Jacobian matrix
J = zeros(2M+1, N+1)
for i=1:2M+1, j=1:N+1
J[i, j] = (i-M-1)^(j-1)
end
e₁ = zeros(N+1)
e₁[1] = 1.0
#Compute filter coefficients
C = J' \ e₁
#Evaluate filter on data matrix
To = typeof(C[1] * one(T)) #Calculate type of output
expr = quote
n = size(data, 1)
smoothed = zeros($To, n)
@inbounds for i in eachindex(smoothed)
smoothed[i] += $(C[M+1])*data[i]
end
smoothed
end
for j=1:M
insert!(expr.args[6].args[3].args[2].args, 1,
:(if i - $j ≥ 1
smoothed[i] += $(C[M+1-j])*data[i-$j]
end)
)
push!(expr.args[6].args[3].args[2].args,
:(if i + $j ≤ n
smoothed[i] += $(C[M+1+j])*data[i+$j]
end)
)
end
return expr
end
It can be used as follows:
window_width = 30
order = 2
svg_filter = SavitzkyGolayFilter{window_with, order}()
filtered_vector= svg_filter(unfiltered_vector)
I have added a line for derivatives computations
struct SavitzkyGolayFilter{M,N,deriv} end
@generated function (::SavitzkyGolayFilter{M,N,deriv})(data::AbstractVector{T}) where {M, N, deriv, T}
#Create Jacobian matrix
J = zeros(2M+1, N+1)
for i=1:2M+1, j=1:N+1
J[i, j] = (i-M-1)^(j-1)/factorial(j-1)
end
e₁ = zeros(N+1)
e₁[deriv+1] = 1.0
#Compute filter coefficients
C = J' \ e₁
#Evaluate filter on data matrix
To = typeof(C[1] * one(T)) #Calculate type of output
expr = quote
n = size(data, 1)
smoothed = zeros($To, n)
@inbounds for i in eachindex(smoothed)
smoothed[i] += $(C[M+1])*data[i]
end
smoothed
end
for j=1:M
insert!(expr.args[6].args[3].args[2].args, 1,
:(if i - $j ≥ 1
smoothed[i] += $(C[M+1-j])*data[i-$j]
end)
)
push!(expr.args[6].args[3].args[2].args,
:(if i + $j ≤ n
smoothed[i] += $(C[M+1+j])*data[i+$j]
end)
)
end
return expr
end
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I was having errors running your code and found this discussion, https://discourse.julialang.org/t/base-call-in-julia-0-6/3983/3 I don't know if you still maintain it but it should be useful for people stumbling on to the same error.