Artifacts: Eye movements, power line noise

Project.toml

@@ -5,13 +5,16 @@ version = "0.4.1"
 
 [deps]
 Artifacts = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
 MixedModels = "ff71e718-51f3-5ec2-a782-8ffcbfa3c316"
 MixedModelsSim = "d5ae56c5-23ca-4a1f-b505-9fc4796fc1fe"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
@@ -23,13 +26,16 @@ ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
 
 [compat]
 Artifacts = "1"
+CSV = "0.10.15"
+CoordinateTransformations = "0.6.4"
 DSP = "0.7,0.8"
 DataFrames = "1"
 Distributions = "0.25"
 FileIO = "1"
 HDF5 = "0.17"
 ImageFiltering = "0.7"
 LinearAlgebra = "1"
+MAT = "0.10.7"
 MixedModels = "4"
 MixedModelsSim = "0.2"
 Parameters = "0.12"

src/AbstractContinuousSignal.jl

@@ -0,0 +1,108 @@
+function simulate_continuoussignal(rng::AbstractRNG, s::EyeMovement, controlsignal::AbstractMatrix, sim::Simulation)
+    # at the lowest level, eyemovement simulation takes a controlsignal = Matrix having size 3 x n_timepoints i.e. a set of gaze direction vectors     
+    headmodel = s.headmodel
+    eye_model = s.eye_model
+    return simulate_eyemovement(headmodel,controlsignal, eye_model)
+end
+
+"""
+controlsignal: actual weights (n_channels x n_timepoints) to be applied to the simulated single-channel powerline noise
+"""
+function simulate_continuoussignal(rng::AbstractRNG, s::PowerLineNoise, controlsignal::AbstractMatrix, sim::Simulation)
+    base_freq = s.base_freq
+    harmonics = s.harmonics
+    sampling_rate = s.sampling_rate
+    weights_harmonics = s.weights_harmonics
+
+    n_samples = size(controlsignal)[end]
+    k = 0:1:n_samples-1
+    # TODO add check for nyquist criterion? -> warn or error?
+
+    harmonics_signals = [sin.(2 * pi * (base_freq.*h)/sampling_rate .* k) for h in harmonics].*weights_harmonics
+
+    return reduce(+,harmonics_signals)' .*controlsignal
+end
+
+
+function simulate_continuoussignal(rng::AbstractRNG, s::AbstractNoise, controlsignal::AbstractArray, sim::Simulation)
+    return reshape(simulate_noise(rng,s,length(controlsignal)),size(controlsignal)) .* controlsignal
+end
+
+function simulate_continuoussignal(rng::AbstractRNG, s::UserDefinedContinuousSignal, controlsignal::AbstractArray, sim::Simulation)
+    return s.signal .* controlsignal # the user has already defined what they want
+end
+
+
+
+# generate_controlsignal - returns controlsignal for the specified type of AbstractContinuousSignal.
+# also takes the simulation object since it might influence the generated controlsignal (e.g. generate blinks right after event A -> controlsignal depends on the design) 
+
+
+# for EyeMovement, always return a 3 x n_timepoints matrix containing gaze direction vectors (n_timepoints number of gaze vectors, in 3 dimensions)
+
+function generate_controlsignal(rng::AbstractRNG, cs::GazeDirectionVectors, sim::Simulation)
+    @assert size(cs.val)[1] == 3 "Please make sure gaze data has the shape 3 x n_timepoints."
+    return cs.val
+end
+
+function generate_controlsignal(rng::AbstractRNG, cs::HREFCoordinates, sim::Simulation)
+    return reduce(hcat,gazevec_from_angle_3d.(cs.val[1,:],cs.val[2,:]))
+end
+
+function generate_controlsignal(rng::AbstractRNG, cs::AbstractContinuousSignal, sim::Simulation)
+    return generate_controlsignal(deepcopy(rng), cs.controlsignal, sim)
+end
+
+ # for PowerLineNoise we do not yet know the required n_timepoints, so we postpone controlsignal generation until the eeg and the known-size artifacts have already been generated
+function generate_controlsignal(rng::AbstractRNG, cs::PowerLineNoise, sim::Simulation)
+    return nothing
+end
+
+# AbstractNoise: returns nothing for now, however in future the controlsignal for AbstractNoise may depend on the Simulation.
+function generate_controlsignal(rng::AbstractRNG, cs::AbstractNoise, sim::Simulation)
+    return nothing
+end
+
+# identity function - in case the user already specified the final controlsignal while creating the artifact
+function generate_controlsignal(rng::AbstractRNG, cs::AbstractMatrix, sim::Simulation)
+    return cs
+end
+
+
+
+# simulate_continououssignal - simulates the continuoussignal based on the signal type and the specified controlsignal.
+
+function simulate_continuoussignal(rng::AbstractRNG, s::Union{ARDriftNoise,DCDriftNoise,LinearDriftNoise}, controlsignal::AbstractMatrix, sim::Simulation)
+     # call the singlechannel simulate_continououssignal(..., controlsignal[i,:]) for each row (channel) of controlsignal
+     return hcat(map(controlsignal_channelwise->simulate_continuoussignal(rng,s,controlsignal_channelwise,sim), eachrow(controlsignal))...)'
+end
+
+function simulate_continuoussignal(rng::AbstractRNG, s::LinearDriftNoise, controlsignal::AbstractVector, sim::Simulation)
+    # simulate single-channel drift
+    n_samples = size(controlsignal)[end]
+    return  range(0,1,length=n_samples).*
+            (rand((rng))*2-1) .*    # slope of the line for this particular channel
+            s.scaling_factor .*     # global scaling factor for the entire current LinearDriftNoise
+            controlsignal           # user-controllable weights (1 x timepoint : current_channel vector taken from the entire artifact controlsignal) 
+end
+
+function simulate_continuoussignal(rng::AbstractRNG, s::ARDriftNoise, controlsignal::AbstractVector, sim::Simulation)
+    n_samples = size(controlsignal)[end]
+    return cumsum(s.σ .* randn((rng),n_samples).*controlsignal) 
+end
+
+function simulate_continuoussignal(rng::AbstractRNG, s::DCDriftNoise, controlsignal::AbstractVector, sim::Simulation)
+    n_samples = size(controlsignal)[end]
+    return ones(n_samples).*rand((rng)).* s.scaling_factor .*controlsignal
+end
+
+# generic drift noise: contains multiple fields having different driftnoise types. Simulate the artifact for each of these.
+function simulate_continuoussignal(rng::AbstractRNG, s::DriftNoise, controlsignal::AbstractMatrix, sim::Simulation)
+    fn = fieldnames(DriftNoise)
+    all_s = []
+    for f in fn
+        push!(all_s,simulate_continuoussignal(rng,getfield(s,f),controlsignal,sim))
+    end
+
+    return sum(all_s)
+end

src/UnfoldSim.jl

@@ -17,6 +17,11 @@ using HDF5, Artifacts, FileIO
 
 using LinearAlgebra # headmodel
 
+using CoordinateTransformations
+import CoordinateTransformations.CartesianFromSpherical # artifacts - gaze direction vector
+using CSV # artifacts - importing sample data from csv file
+using MAT # artifacts - loading model from .mat file
+
 import DSP.hanning
 import Base.length
 import Base.size
@@ -31,6 +36,8 @@ include("predefinedSimulations.jl")
 include("headmodel.jl")
 include("helper.jl")
 include("bases.jl")
+include("AbstractContinuousSignal.jl")
+include("UnfoldSimArtifacts.jl")
 
 export size, length
 export AbstractComponent, AbstractNoise, AbstractOnset, AbstractDesign
@@ -78,3 +85,10 @@ export AbstractHeadmodel, Hartmut, headmodel, leadfield, orientation, magnitude
 # multichannel
 export MultichannelComponent
 end
+
+# artifacts
+export simulate_eyemovement 
+export AbstractControlSignal, GazeDirectionVectors, HREFCoordinates
+export AbstractContinuousSignal, EyeMovement, PowerLineNoise
+export az_simulation
+export example_data_eyemovements