diff --git a/Project.toml b/Project.toml
index f734b2a..b0c7639 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "KissSmoothing"
 uuid = "23b0397c-cd08-4270-956a-157331f0528f"
 authors = ["Francesco Alemanno <francescoalemanno710@gmail.com>"]
-version = "1.0.6"
+version = "1.0.7"
 
 [deps]
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
diff --git a/README.md b/README.md
index fd423c3..5ed7179 100644
--- a/README.md
+++ b/README.md
@@ -154,15 +154,19 @@ savefig("nspline.png")
 
     fit_sine_series(X::Vector, Y::Vector, basis_elements::Integer, noise=0)
 
-fit Y ~ 1 + X + Σ sin(.) according to:
+fit Y ~ 1 + X + Σ sin(.) by minimising Σ (Y - f(x))^2 + lambda * ∫(Δ^order F)^2
 
-    `X` : array N, N number of training points
+    `X` : array N, N number of training points.
 
-    `Y` : array N, N number of training points
+    `Y` : array N, N number of training points.
 
-    `basis_elements` : number of sine terms
+    `basis_elements` : number of sine terms.
     
-    `noise` : noise filtering according to Wiener method, default to 0 (off)
+    Keyword arguments:
+    
+    `lambda` : intensity of regularization.
+    
+    `order` : derivative to regularise.
 
 returns a callable function.
 
@@ -170,17 +174,20 @@ returns a callable function.
 
 ```julia
 using PyPlot, KissSmoothing
-t = LinRange(0,pi,250)
-ty = sin.(t.^2) .+ t
-y = ty .+ randn(length(t)) .*0.05
-fn = fit_sine_series(t,y,20)
-scatter(t, y, color="gray",s=2,label="noisy")
-plot(t, fn.(t), color="red",lw=1.,label="sine fit")
-plot(t,ty, color="blue",lw=0.7,label="true")
+fg(x) = sin(x^2) + x
+x = collect(LinRange(0,pi,250))
+xc = identity.(x)
+filter!(x->(x-1)^2>0.1, x)
+filter!(x->(x-2)^2>0.1, x)
+y = fg.(x) .+ randn(length(x)) .*0.05
+fn = fit_sine_seris(x,y,20, lambda = 0.00001)
+scatter(x, y, color="gray",s=2,label="noisy")
+plot(xc, fn.(xc), color="red",lw=1.,label="fit")
+plot(xc,fg.(xc), color="blue",lw=0.7,label="true")
 xlabel("X")
 ylabel("Y")
 legend()
 tight_layout()
 savefig("sine_fit.png")
 ```
-![sine_fit.png](sine_fit.png "Plot of NSpline estimation")
+![sine_fit.png](sine_fit.png "Plot of FFT Spline estimation")
diff --git a/sine_fit.png b/sine_fit.png
index d8cf120..384320a 100644
Binary files a/sine_fit.png and b/sine_fit.png differ
diff --git a/src/KissSmoothing.jl b/src/KissSmoothing.jl
index 7052e66..ecd2392 100644
--- a/src/KissSmoothing.jl
+++ b/src/KissSmoothing.jl
@@ -198,19 +198,23 @@ end
 """
     fit_sine_series(X::Vector, Y::Vector, basis_elements::Integer, noise=0)
 
-fit Y ~ 1 + X + Σ sin(.) according to:
+fit Y ~ 1 + X + Σ sin(.) by minimising Σ (Y - f(x))^2 + lambda * ∫(Δ^order F)^2
 
-    `X` : array N, N number of training points
+    `X` : array N, N number of training points.
 
-    `Y` : array N, N number of training points
+    `Y` : array N, N number of training points.
 
-    `basis_elements` : number of sine terms
+    `basis_elements` : number of sine terms.
     
-    `noise` : noise filtering according to Wiener method, default to 0 (off)
+    Keyword arguments:
+    
+    `lambda` : intensity of regularization.
+    
+    `order` : derivative to regularise.
 
 returns a callable function.
 """
-function fit_sine_series(X::AbstractVector{<:Real}, Y::AbstractVector{<:Real}, basis_elements::Integer, noise = 0.0)
+function fit_sine_series(X::AbstractVector{<:Real},Y::AbstractVector{<:Real}, basis_elements::Integer; lambda = 0.0, order::Int64=3)
     lx, hx = extrema(X)
     T = @. (X - lx)/(hx-lx)*pi
     M = zeros(length(X),2+basis_elements)
@@ -218,22 +222,24 @@ function fit_sine_series(X::AbstractVector{<:Real}, Y::AbstractVector{<:Real}, b
         M[i, 1] = 1
         M[i, 2] = T[i]
         for k in 1:basis_elements
-            M[ i, 2+k] = sin(k*T[i])
+            M[i, 2+k] = sin(k*T[i])#/k^order
         end
     end
-    C = M\Y
+    phi = M'M
+    for i=3:2+basis_elements
+        phi[i,i] += lambda*(i-2)^(2order)
+    end
+    C = phi \ (M'*Y)
     return function fn(x)
         t = (x - lx)/(hx-lx)*pi
-        s = C[1]+C[2]*t
+        s = C[1] + C[2]*t
         for k in 1:basis_elements
-            cn1 = C[2+k]
-            SN = abs2(cn1)
-            hn = max(1 - noise*noise/SN,0)
-            s += cn1*sin(k*t)*hn
+            s += C[2+k]*sin(k*t)#/k^order
         end
         s
     end
 end
+
     
 export denoise, fit_rbf, RBF, fit_nspline, fit_sine_series
 end # module
diff --git a/test/runtests.jl b/test/runtests.jl
index 86a9ad5..4368ba1 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -60,7 +60,7 @@ end
     for μ in LinRange(-100,100,5)
         t = LinRange(0,2pi,150)
         y = sin.(t) .+ μ .* t
-        fn = fit_sine_series(t,y,50)
+        fn = fit_sine_series(t,y,50, order = 3, lambda = 0.00001)
         pred_y = fn.(t)
         error = sqrt(sum(abs2, pred_y .- y)/length(t))
         @test error < 0.0002