Skip to content

Instantly share code, notes, and snippets.

@jiahao
Last active April 26, 2024 12:03
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 2 You must be signed in to fork a gist
  • Save jiahao/b8b5ac328c18b7ae8a17 to your computer and use it in GitHub Desktop.
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
"""
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;])
@j605
Copy link

j605 commented Oct 26, 2017

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.

@Helveg
Copy link

Helveg commented Jun 9, 2019

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)

@mrozkamil
Copy link

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