FIR with duration

src/basisfunctions.jl

@@ -83,9 +83,14 @@ end
 $(SIGNATURES)
 Generate a sparse FIR basis around the *τ* timevector at sampling rate *sfreq*. This is useful if you cannot make any assumptions on the shape of the event responses. If unrounded events are supplied, they are split between samples. E.g. event-latency = 1.2 will result in a "0.8" and a "0.2" entry.
 
+
+Advanced: second input can be duration in samples - careful: `times(firbasis)` always assumes duration = 1. Therefore, 
+issues with LMM and predict will appear!
+
 # keyword arguments
 `interpolate` (Bool, default false): if true, interpolates events between samples linearly. This results in `predict` functions to return a trailling 0
 
+
 # Examples
 Generate a FIR basis function from -0.1s to 0.3s at 100Hz
 ```julia-repl
@@ -97,12 +102,15 @@ julia>  f(103.3)
 ```
 
 """
-function firbasis(τ, sfreq, name = ""; interpolate = false)
+function firbasis(τ, sfreq, name = ""; interpolate = false, max_height=nothing)
     τ = round_times(τ, sfreq)
     if interpolate
         # stop + 1 step, because we support fractional event-timings
         τ = (τ[1], τ[2] + 1 ./ sfreq)
     end
+  if !isnothing(max_height)
+        τ = (τ[1], τ[2] + max_height)
+  end
     times = range(τ[1], stop = τ[2], step = 1 ./ sfreq)
 
     shift_onset = Int64(floor(τ[1] * sfreq))
@@ -119,34 +127,57 @@ firbasis(; τ, sfreq, name = "", kwargs...) = firbasis(τ, sfreq, name; kwargs..
 """
 $(SIGNATURES)
 Calculate a sparse firbasis
+
+second input can be duration in samples - careful: `times(firbasis)` always assumes duration = 1. Therefore, 
+issues with LMM and predict will appear!
+
 # Examples
 
 ```julia-repl
 julia>  f = firkernel(103.3,range(-0.1,step=0.01,stop=0.31))
+julia>  f_dur = firkernel([103.3 4],range(-0.1,step=0.01,stop=0.31))
 ```
 """
-function firkernel(e, times; interpolate = false)
-    @assert ndims(e) <= 1 #either single onset or a row vector where we will take the first one :)
-    if size(e, 1) > 1
-        # XXX we will soon assume that the second entry would be the duration
-        e = Float64(e[1])
-    end
-    ksize = length(times) # kernelsize
-    if interpolate
-        eboth = [1 .- (e .% 1) e .% 1]
-        eboth[isapprox.(eboth, 0, atol = 1e-15)] .= 0
-        return spdiagm(
-            ksize + 1,
+function firkernel(ev, times; interpolate = false)
+    @assert ndims(ev) <= 1 #either single onset or a row vector where we will take the first one :)
+
+    # duration is 1 for FIR and duration is duration (in samples!) if e is vector
+    e = ev[1] # latency in samples  
+    dur = size(ev, 1) == 1 ? 1 : Int.(ceil(ev[2]))
+
+    ksize = length(times) #isnothing(maxsize) ? length(times) : max_height # kernelsize
+
+    # if interpolatethe first and last entry is split, depending on whether we have "half" events
+    eboth = interpolate ? [1 .- (e .% 1) e .% 1] : [e]
+
+     values =[eboth[1] ; repeat([1],dur-1); eboth[2:end]] 
+     values[isapprox.(eboth, 0, atol = 1e-15)] .= 0 # keep sparsity pattern
+
+     # without interpolation, remove last entry again, which should be 0 anyway
+  # values = interpolate ? values : values[1:end-1]  # commented out if the eboth[2:end] trick works
+
+        # build the matrix, we define it by diagonals which go from 0, -1, -2 ...
+    pairs = [x => repeat([y],ksize) for (x,y) = zip(.-range(0,dur),values)]
+            return spdiagm(
+            ksize + interpolate ? 1 : 0,
             ksize,
-            0 => repeat([eboth[1]], ksize),
-            -1 => repeat([eboth[2]], ksize),
+          pairs...
         )
-    else
-        #eboth = Int(round(e))
-        return spdiagm(ksize, ksize, 0 => repeat([1], ksize))
+
     end
 
 
+splinebasis(; τ, sfreq, nsplines, name) = splinebasis(τ, sfreq, nsplines, name)
+splinebasis(τ, sfreq, nsplines) = splinebasis(τ, sfreq, nsplines, "basis_" * string(rand(1:10000)))
+
+shiftOnset(basis::Union{SplineBasis,FIRBasis}) = basis.shiftOnset
+
+function splinebasis(τ, sfreq, nsplines, name::String)
+    τ = Unfold.round_times(τ, sfreq)
+    times = range(τ[1], stop = τ[2], step = 1 ./ sfreq)
+    kernel = e -> splinekernel(e, times, nsplines - 2)
+
+
 end
 
 

test/basisfunctions.jl

@@ -24,6 +24,14 @@
     @test Unfold.kernel(firbase_on, 0.5)[1:3, 1:3] ==
           [0.5 0.0 0.0; 0.5 0.5 0.0; 0.0 0.5 0.5]
 
+# defining a firkernel with a duration = 1 sample, should be equal to the non-duration one
+f_dur = Unfold.firkernel([103.3;1],range(-0.1,step=0.01,stop=0.31))
+f_fir = Unfold.firkernel(103.3,range(-0.1,step=0.01,stop=0.31))
+@test f_dur == f_fir
+
+# test duration for samples = 4
+f_dur = Unfold.kernel(firbase)([1,4])
+@test all(sum(f_dur,dims=1) .== 4)
 end
 @testset "BOLD" begin