{{ message }}

Instantly share code, notes, and snippets.

# jiahao/savitzkygolay.jl

Last active Jul 20, 2021
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.0 #Compute filter coefficients C = J' \ e₁ #Evaluate filter on data matrix To = typeof(C * 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.args.args.args, 1, :(if i - \$j ≥ 1 smoothed[i] += \$(C[M+1-j])*data[i-\$j] end) ) push!(expr.args.args.args.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 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 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.0 #Compute filter coefficients C = J' \ e₁ #Evaluate filter on data matrix To = typeof(C * 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.args.args.args, 1, :(if i - \$j ≥ 1 smoothed[i] += \$(C[M+1-j])*data[i-\$j] end) ) push!(expr.args.args.args.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) ``````