Блог

Я наткнулся на реализацию SG-filter в Julia по этой ссылке. Когда я выполняю функцию apply_filter , возвращается ошибка —

 UndefVarError: apply_filter not defined
  

Я думаю, что это реализация для предыдущей версии Julia (?). На данный момент я выполняю это в Julia 1.0. Не удалось найти документацию об определенных типах, в чем, как я предполагаю, заключается ошибка

Я хотел бы предупредить пользователя об использовании функции savitzkyGolay в Julia. Есть несоответствие с результатом реализации Scipy (который, должно быть, прошел несколько итераций проверки сообществом)

 @pyimport scipy.signal as ss
x=[1,2,3,4,5,6,7,8,9,10]


savitzkyGolay(x,5,1)

10-element Array{Float64,1}:
 1.6000000000000003
 2.200000000000001 
 3.0               
 4.0               
 5.000000000000001 
 6.000000000000001 
 7.0               
 8.0               
 8.8               
 9.400000000000002 


#Python's scipy implementation
ss.savgol_filter(x,5,1)

10-element Array{Float64,1}:
 1.0000000000000007
 2.0000000000000004
 2.9999999999999996
 3.999999999999999 
 4.999999999999999 
 5.999999999999999 
 6.999999999999998 
 7.999999999999998 
 8.999999999999996 
 9.999999999999995 
  

Если это может помочь, я упростил код.

 using Pkg, LinearAlgebra, DSP, Plots
function vandermonde(halfWindow, polyDeg)
    x=[1.0*i for i in -halfWindow:halfWindow]

    n = polyDeg 1
    m = length(x)

    V = zeros(m, n)

    for i = 1:m
        V[i,1] = 1.0
    end
    for j = 2:n
        for i = 1:m
            V[i,j] = x[i] * V[i,j-1]
        end
    end

    return V
end

function SG(halfWindow, polyDeg)

    V = vandermonde(halfWindow,polyDeg)
    Q,R=qr(V)
    n = polyDeg 1
    m = 2*halfWindow 1
    R1 = vcat(R, zeros(m-n,n))
    sg = R1Q'

    for i in 1:(polyDeg 1)

        sg[i,:] = sg[i,:]*factorial(i-1)

    end

    return sg'
end

function apply_filter(filter,signal)

    halfWindow = round(Int,(length(filter)-1)/2)

    padded_signal = [signal[1]*ones(halfWindow);signal;signal[end]*ones(halfWindow)]

    filter_cross_signal = conv(filter[end:-1:1], padded_signal)

    return filter_cross_signal[2*halfWindow 1:end-2*halfWindow]
end
  

Вот как я это использую :

 mean_speed_unfiltered = readdlm("mean_speeds_raw_-2.txt")
sg = SG(500,2); # halt-window, polynomal degree

t =  10*10^(-3)#s #time of the simulation
dt = 0.1/γ; #time step
Nt = convert(Int, round(t/dt)); #number of iteration

#Smooth the mean speed curve: 

mean_speeds_smoothed = apply_filter(sg[:,1],mean_speed_unfiltered)

png(plot([j*dt for j=0:Nt] , mean_speeds_smoothed, title = "Smoothed mean speed over 
time", xlabel = "t (s)"), "Mean_speed_filtered_SG")

derivative_mean_speeds_smoothed = apply_filter(sg[:,2],mean_speed_unfiltered)

plt1 = plot(mean_speeds_smoothed,derivative_mean_speeds_smoothed, title = "derivative mean speed over speed", xlabel = "<v>(t) (s)", ylabel = "d<v(t)>/dt")
png(plt1, "Force_SG_1D2Lasers")
  

Однако мне кажется, что код, представленный в https://gist.github.com/lnacquaroli/c97fbc9a15488607e236b3472bcdf097#file-savitzkygolay-jl-L34 быстрее.