add fit_sine_series

* add fit_sine_series, with wiener filtering

Project.toml

@@ -1,7 +1,7 @@
 name = "KissSmoothing"
 uuid = "23b0397c-cd08-4270-956a-157331f0528f"
 authors = ["Francesco Alemanno <[email protected]>"]
-version = "1.0.5"
+version = "1.0.6"
 
 [deps]
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

README.md

@@ -149,3 +149,38 @@ tight_layout()
 savefig("nspline.png")
 ```
 ![nspline.png](nspline.png "Plot of NSpline estimation")
+
+## Sine Series Estimation
+
+    fit_sine_series(X::Vector, Y::Vector, basis_elements::Integer, noise=0)
+
+fit Y ~ 1 + X + Σ sin(.) according to:
+
+    `X` : array N, N number of training points
+
+    `Y` : array N, N number of training points
+
+    `basis_elements` : number of sine terms
+
+    `noise` : noise filtering according to Wiener method, default to 0 (off)
+
+returns a callable function.
+
+### Example
+
+```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")
+xlabel("X")
+ylabel("Y")
+legend()
+tight_layout()
+savefig("sine_fit.png")
+```
+![sine_fit.png](sine_fit.png "Plot of NSpline estimation")

src/KissSmoothing.jl

@@ -195,6 +195,45 @@ function fit_nspline(
     end
 end
 
+"""
+    fit_sine_series(X::Vector, Y::Vector, basis_elements::Integer, noise=0)
+
+fit Y ~ 1 + X + Σ sin(.) according to:
+
+    `X` : array N, N number of training points
+
+    `Y` : array N, N number of training points
 
-export denoise, fit_rbf, RBF, fit_nspline
+    `basis_elements` : number of sine terms
+    
+    `noise` : noise filtering according to Wiener method, default to 0 (off)
+
+returns a callable function.
+"""
+function fit_sine_series(X::AbstractVector{<:Real}, Y::AbstractVector{<:Real}, basis_elements::Integer, noise = 0.0)
+    lx, hx = extrema(X)
+    T = @. (X - lx)/(hx-lx)*pi
+    M = zeros(length(X),2+basis_elements)
+    for i in eachindex(X)
+        M[i, 1] = 1
+        M[i, 2] = T[i]
+        for k in 1:basis_elements
+            M[ i, 2+k] = sin(k*T[i])
+        end
+    end
+    C = M\Y
+    return function fn(x)
+        t = (x - lx)/(hx-lx)*pi
+        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
+        end
+        s
+    end
+end
+
+export denoise, fit_rbf, RBF, fit_nspline, fit_sine_series
 end # module

test/runtests.jl

@@ -55,6 +55,18 @@ end
     end
 end
 
+
+@testset "Fit Sine Series" begin
+    for μ in LinRange(-100,100,5)
+        t = LinRange(0,2pi,150)
+        y = sin.(t) .+ μ .* t
+        fn = fit_sine_series(t,y,50)
+        pred_y = fn.(t)
+        error = sqrt(sum(abs2, pred_y .- y)/length(t))
+        @test error < 0.0002
+    end
+end
+
 @testset "Fit NSpline" begin
     for μ in LinRange(-100,100,5)
         t = LinRange(0,2pi,150)