@nsajko Thank you again for your help. This is my current implementation of Savitzky-Golay filter:

```
struct SavitzkyGolayFilterK{T1 <: Int, T2 <: AbstractFloat}
window_size :: T1 # window size must be an odd number
M :: T1 # window half-size must be an odd number
degree :: T1 # degree of the fitting polynomial
deriv :: T1 # order of the derivative to compute
dx::T2 # grid spacing
C :: Vector{T2} # filter coefficients
proj_matrix:: Matrix{T2} # projection matrix
function SavitzkyGolayFilterK(window_size::T1, degree::T1, deriv::T1, dx::T2 ) where {T1, T2}
isodd(window_size) || throw(ArgumentError("window_size must be an odd number."))
window_size >= degree + 1 || throw(ArgumentError("window_size must be >= degree + 1."))
M = (window_size - 1) ÷ 2
#Create Jacobian matrix
J = zeros(T2,(2M+1, degree+1))
for i=1:2M+1, j=1:degree+1
J[i, j] = (i-M-1)^(j-1)
end
proj_mat = (J' * J) \ J'
C = proj_mat[deriv+1,:] .* T2(factorial(deriv+1))
C = C ./ dx^deriv
return new{T1,T2}(window_size, M, degree, deriv, dx, C, proj_mat)
end
end
SavitzkyGolayFilterK(window_size, degree) = SavitzkyGolayFilterK(window_size, degree, 0, 1.0)
SavitzkyGolayFilterK(window_size, degree, deriv) = SavitzkyGolayFilterK(window_size, degree, deriv, 1.0)
function (sg_filter::SavitzkyGolayFilterK)(data::Vector{T}; mode = :wrap) where T
M = sg_filter.M
if T != typeof(sg_filter.C[1])
C = T.(sg_filter.C)
else
C = @view sg_filter.C[:]
end
n = size(data, 1)
smoothed = C[M+1] .* data
if mode == :wrap
@inbounds for i in eachindex(smoothed)
for j=1:M
smoothed[i] += (C[M+1-j]) * (i - j ≥ 1 ? data[i-j] : data[i-j+n])
smoothed[i] += (C[M+1+j]) * (i + j ≤ n ? data[i+j] : data[i+j-n])
end
end
elseif mode == :reflect
@inbounds for i in eachindex(smoothed)
for j=1:M
smoothed[i] += (C[M+1-j]) * (i - j ≥ 1 ? data[i-j] : data[1-(i-j)])
smoothed[i] += (C[M+1+j]) * (i + j ≤ n ? data[i+j] : data[n-(i+j)])
end
end
elseif mode == :constant
@inbounds for i in eachindex(smoothed)
for j=1:M
smoothed[i] += (C[M+1-j]) * (i - j ≥ 1 ? data[i-j] : data[1])
smoothed[i] += (C[M+1+j]) * (i + j ≤ n ? data[i+j] : data[n])
end
end
elseif mode == :pad
@inbounds for i in eachindex(smoothed)
for j=1:M
smoothed[i] += (C[M+1-j]) * (i - j ≥ 1 ? data[i-j] : 0)
smoothed[i] += (C[M+1+j]) * (i + j ≤ n ? data[i+j] : 0)
end
end
end
return smoothed
end
```

I call it by defining the structure first, with e.g.

```
signal = rand(1024)
window_length = 17
M = 8
N = 3
deriv = 1
SG_filterK = SavitzkyGolayFilterK(window_length, N, deriv)
```

I benchmarked my code against other implementations

```
Savitzky-Golay filter form SavitzkyGolay.jl
21.558 μs (66 allocations: 112.09 KiB)
Savitzky-Golay filter form SavitzkyGolay.jl with a constructor builded beforehand
21.101 μs (63 allocations: 103.94 KiB)
my Savitzky-Golay filter
16.775 μs (1 allocation: 8.12 KiB)
Savitzky-Golay filter of jiahao
9.506 μs (1 allocation: 8.12 KiB)
```

I’m a bit disappointed that my implementation is 2 times slower than jiahao’s, however I managed to modify his code to behave better at the edges

```
struct SavitzkyGolayFilter{M,N,deriv, mode} end
@generated function (::SavitzkyGolayFilter{M,N,deriv,mode})(data::Vector{T}) where {M, N, deriv, mode, T}
#Create Jacobian matrix
J = zeros(T,(2M+1, N+1))
for i=1:2M+1, j=1:N+1
J[i, j] = (i-M-1)^(j-1)
end
proj_mat = (J' * J) \ J'
C = proj_mat[deriv+1,:] .* T(factorial(deriv+1))
#Evaluate filter on data matrix
expr = quote
n = size(data, 1)
# smoothed = zeros($To, n)
smoothed = zeros(T, n)
@inbounds for i in eachindex(smoothed)
smoothed[i] += $(C[M+1])*data[i]
end
smoothed
end
for j=1:M
if mode == :constant
tmp_expr = :(if i - $j ≥ 1
smoothed[i] += $(C[M+1-j])*data[i-$j]
else
smoothed[i] += $(C[M+1-j])*data[1]
end)
elseif mode == :wrap
tmp_expr = :(if i - $j ≥ 1
smoothed[i] += $(C[M+1-j])*data[i-$j]
else
smoothed[i] += $(C[M+1-j])*data[i-$j+n]
end)
elseif mode == :reflect
tmp_expr = :(if i - $j ≥ 1
smoothed[i] += $(C[M+1-j])*data[i-$j]
else
smoothed[i] += $(C[M+1-j])*data[1-(i-$j)]
end)
elseif mode == :pad
tmp_expr = :(if i - $j ≥ 1
smoothed[i] += $(C[M+1-j])*data[i-$j]
end)
end
insert!(expr.args[6].args[3].args[2].args, 1, tmp_expr)
if mode == :constant
tmp_expr = :(if i + $j ≤ n
smoothed[i] += $(C[M+1+j])*data[i+$j]
else
smoothed[i] += $(C[M+1+j])*data[n]
end)
elseif mode == :reflect
tmp_expr = :(if i + $j ≤ n
smoothed[i] += $(C[M+1+j])*data[i+$j]
else
smoothed[i] += $(C[M+1+j])*data[2*n-(i+$j)]
end)
elseif mode == :wrap
tmp_expr = :(if i + $j ≤ n
smoothed[i] += $(C[M+1+j])*data[i+$j]
else
smoothed[i] += $(C[M+1+j])*data[i+$j-n]
end)
else
tmp_expr = :(if i + $j ≤ n
smoothed[i] += $(C[M+1+j])*data[i+$j]
end)
end
push!(expr.args[6].args[3].args[2].args,tmp_expr)
end
return expr
end
```