ICA deconstruct

NeuroAnalyzer tutorials: Process EEG (5)

Load data:

using NeuroAnalyzer
using Plots
eeg = load("files/eeg.hdf");
trim!(eeg, seg=(0, 20), remove_epochs=false)

[ Info: Precompiling NeuroAnalyzer [b40dcafa-b65c-47ad-a230-3174ba5cadd1] (cache misses: include_dependency fsize change (2), wrong dep version loaded (2), incompatible header (4), mismatched flags (10))

[ Info: NeuroAnalyzer v0.25.5-dev
[ Info: NeuroAnalyzer path: /home/eb/Documents/Code/NeuroAnalyzer.jl
[ Info:  Preferences:
[ Info:     Use CUDA: false
[ Info: Progress bar: true
[ Info:      Verbose: true
[ Info: Exclude bads: false
[ Info: Preparing resources
[ Info: Loading plugins:
[ Info:  Loaded: na_test_plugin.jl
[ Info:  Loaded: plot_env.jl
[ Info:  Loaded: plot_ispc.jl
[ Info:  Loaded: plot_itpc.jl
[ Info:  Loaded: plot_pli.jl
[ Info: Loaded: EEG (24 × 282991 × 1; 1105.43 s)

Deconstruct all EEG channels into defined number (3) of components:

ic, ic_mw, ic_var = ica_decompose(eeg, ch="eeg", n=3)

[ Warning: The input signal should be cleaned from major artifacts and HP filtered at 1-2 Hz prior to ICA decomposition.
[ Info: Attempting to calculate 3 components
[ Info: Training will end when W change = 0.99 or after 900 steps
[ Info: Data will be demeaned and pre-whitened
[ Info: Converged at: 1.0e-6
[ Info: Component  1: percent variance accounted for: 54.51
[ Info: Component  2: percent variance accounted for: 25.37
[ Info: Component  3: percent variance accounted for: 6.39

(ic = [0.041414342896599704 0.20233594094662544 … 0.1307358949251824 0.003508170307992964; -0.03447907823662653 -0.08877993613537007 … -0.043490062850950736 -0.0011364618553820073; 0.2509040892011964 0.20955126613223773 … -0.056342466250726544 -0.002197029501350737],
 ic_mw = [-0.7291382188128523 1.1840190630854117 -10.488505796390916; -0.717872833610213 -2.408176101773287 -8.464478486652192; … ; -0.04264383161788086 -1.0608777101229097 -9.281003690985058; 0.09672739692726243 1.3238328346662038 -13.522009529774124],
 ic_var = [54.513912049710754, 25.374017851427155, 6.3933170467389004],)

(!) Components are sorted in the order of decreasing variance accounted for (the higher the value, the higher the component accounts for all the data).

(!) By default, 100 iterations per each tolerance value ([0.000001, 0.00001, 0.0001, 0.001, 0.01, 0.1, 0.5, 0.9, 0.99]); hence the default 100 iterations gives 900 steps.

ICA components and weights may be included in the NeroAnalyzer object as embedded components:

eeg_ic = deepcopy(eeg)
add_component!(eeg_ic, c=:ic, v=ic)
add_component!(eeg_ic, c=:ic_mw, v=ic_mw)

ICA reconstruct

Reconstruct the source signal without a specific component(s):

ica_reconstruct(eeg, ic, ic_mw, ch="eeg", ic_idx=[1, 3])

(!) By default ic_idx is the list of components that should be removed during reconstruction.

Reconstruct the source signal keeping only the selected component(s), use keep=true option:

ica_reconstruct(eeg, ic, ic_mw, ch="eeg", ic_idx=1:2, keep=true)
# is equivalent to
ica_reconstruct(eeg, ic, ic_mw, ch="eeg", ic_idx=3, keep=false)

Use the embedded component for reconstruction:

ica_reconstruct(eeg_ic, ch="eeg", ic_idx=1)

NeuroAnalyzer.NEURO(NeuroAnalyzer.HEADER(Dict{Symbol, Any}(:weight => -1, :id => "", :middle_name => "", :height => -1, :head_circumference => -1, :handedness => "", :last_name => "528004  SIT 52, 20220831-122227-{d589f756-53fc-4f1b-915d-6e3b8c1560ad}", :first_name => ""), Dict{Symbol, Any}(:epoch_id => "", :channel_type => ["eeg", "eeg", "eeg", "eeg", "eeg", "eeg", "eeg", "eeg", "eeg", "eeg"  …  "eeg", "eeg", "ref", "ref", "eeg", "eeg", "eeg", "eog", "eog", "ecg"], :label => ["Fp1", "Fp2", "F3", "F4", "C3", "C4", "P3", "P4", "O1", "O2"  …  "T5", "T6", "A1", "A2", "Fz", "Cz", "Pz", "EOG1", "EOG2", "ECG"], :prefiltering => ["HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz"  …  "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz"], :gain => [0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218  …  0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218], :data_type => "eeg", :recording_notes => "", :recording_date => "31.08.22", :sampling_rate => 256, :file_type => "EDF"…), Dict(:name => "", :design => "", :notes => "")), [0.0, 0.004, 0.008, 0.012, 0.016, 0.02, 0.023, 0.027, 0.031, 0.035  …  1085.391, 1085.395, 1085.398, 1085.402, 1085.406, 1085.41, 1085.414, 1085.418, 1085.422, 1085.426], [0.0, 0.004, 0.008, 0.012, 0.016, 0.02, 0.023, 0.027, 0.031, 0.035  …  1085.391, 1085.395, 1085.398, 1085.402, 1085.406, 1085.41, 1085.414, 1085.418, 1085.422, 1085.426], [-2.672432879834711 -2.3029968062728146 … 0.5394552203833965 0.021697964158517335; -2.040740573035974 -1.5599444635088975 … 0.5816413234864798 0.021333509229431388; … ; 6.441793000605642 5.047332343535441 … -0.6385598239986807 0.0008236881401195717; -27.729377847865074 -33.98805752002153 … -8.80813398237423 0.0018241535887852933;;;], Dict{Any, Any}(), 0×5 DataFrame
 Row │ id      start    length   description  channel 
     │ String  Float64  Float64  String       Int64   
─────┴────────────────────────────────────────────────, 23×9 DataFrame
 Row │ label   loc_radius  loc_theta  loc_x    loc_y    loc_z    loc_radius_sp ⋯
     │ String  Float64     Float64    Float64  Float64  Float64  Float64       ⋯
─────┼──────────────────────────────────────────────────────────────────────────
   1 │ Fp1           1.0       108.0    -0.31     0.95    -0.03            1.0 ⋯
   2 │ Fp2           1.0        72.0     0.31     0.95    -0.03            1.0
   3 │ F7            1.0       144.0    -0.81     0.59    -0.03            1.0
   4 │ F3            0.65      129.0    -0.55     0.67     0.5             1.0
   5 │ Fz            0.51       90.0     0.0      0.72     0.7             1.0 ⋯
   6 │ F4            0.65       51.0     0.55     0.67     0.5             1.0
   7 │ F8            1.0        36.0     0.81     0.59    -0.03            1.0
   8 │ T3            1.0       180.0    -1.0      0.0     -0.03            1.0
   9 │ C3            0.51      180.0    -0.72     0.0      0.7             1.0 ⋯
  10 │ Cz            0.0         0.0     0.0      0.0      1.0             1.0
  11 │ C4            0.51        0.0     0.72     0.0      0.7             1.0
  ⋮  │   ⋮         ⋮           ⋮         ⋮        ⋮        ⋮           ⋮       ⋱
  14 │ P3            0.65      231.0    -0.55    -0.67     0.5             1.0
  15 │ Pz            0.51      270.0     0.0     -0.72     0.7             1.0 ⋯
  16 │ P4            0.65      309.0     0.55    -0.67     0.5             1.0
  17 │ T6            1.0       324.0     0.81    -0.59    -0.03            1.0
  18 │ O1            1.0       252.0    -0.31    -0.95    -0.03            1.0
  19 │ O2            1.0       288.0     0.31    -0.95    -0.03            1.0 ⋯
  20 │ A1            1.0       192.0    -0.92    -0.23    -0.55            1.1
  21 │ A2            1.0       -12.0     0.92    -0.23    -0.55            1.1
  22 │ EOG1          1.01      149.0    -0.87     0.51    -0.37            1.0
  23 │ EOG2          1.01       31.0     0.87     0.51    -0.37            1.0 ⋯
                                                    3 columns and 2 rows omitted, ["reset_components(OBJ)", "filter(OBJ, ch=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24], fprototype=iirnotch, ftype=nothing, cutoff=50, order=8, rp=-1, rs=-1, dir=twopass, w=nothing)", "reset_components(OBJ)", "filter(OBJ, ch=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24], fprototype=fir, ftype=hp, cutoff=0.1, order=8, rp=-1, rs=-1, dir=twopass, w=nothing)", "reset_components(OBJ)", "trim(OBJ, seg=(1, 257), remove_epochs=false)", "reset_components(OBJ)", "trim(OBJ, seg=(1793, 2561), remove_epochs=false)", "reset_components(OBJ)", "trim(OBJ, seg=(1537, 2049), remove_epochs=false)"  …  "reset_components(OBJ)", "trim(OBJ, seg=(17409, 19585), remove_epochs=false)", "reset_components(OBJ)", "trim(OBJ, seg=(49537, 51265), remove_epochs=false)", "reset_components(OBJ)", "trim(OBJ, seg=(1, 5121), remove_epochs=false)", "add_component(OBJ, c=ic, v", "add_component(OBJ, c=ic_mw, v", "reset_components(OBJ)", "ica_reconstruct(OBJ, ch=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 19, 20, 21], ic_idx=1, keep=false)"])

(!) The components must be named :ic and :ic_mw.

Reconstruct using external components:

ica_remove(eeg, ic, ic_mw, ch="eeg", ic_idx=1)

NeuroAnalyzer.NEURO(NeuroAnalyzer.HEADER(Dict{Symbol, Any}(:weight => -1, :id => "", :middle_name => "", :height => -1, :head_circumference => -1, :handedness => "", :last_name => "528004  SIT 52, 20220831-122227-{d589f756-53fc-4f1b-915d-6e3b8c1560ad}", :first_name => ""), Dict{Symbol, Any}(:epoch_id => "", :channel_type => ["eeg", "eeg", "eeg", "eeg", "eeg", "eeg", "eeg", "eeg", "eeg", "eeg"  …  "eeg", "eeg", "ref", "ref", "eeg", "eeg", "eeg", "eog", "eog", "ecg"], :label => ["Fp1", "Fp2", "F3", "F4", "C3", "C4", "P3", "P4", "O1", "O2"  …  "T5", "T6", "A1", "A2", "Fz", "Cz", "Pz", "EOG1", "EOG2", "ECG"], :prefiltering => ["HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz"  …  "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz", "HP:0,18Hz LP:104,0Hz"], :gain => [0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218  …  0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218, 0.17935713740749218], :data_type => "eeg", :recording_notes => "", :recording_date => "31.08.22", :sampling_rate => 256, :file_type => "EDF"…), Dict(:name => "", :design => "", :notes => "")), [0.0, 0.004, 0.008, 0.012, 0.016, 0.02, 0.023, 0.027, 0.031, 0.035  …  1085.391, 1085.395, 1085.398, 1085.402, 1085.406, 1085.41, 1085.414, 1085.418, 1085.422, 1085.426], [0.0, 0.004, 0.008, 0.012, 0.016, 0.02, 0.023, 0.027, 0.031, 0.035  …  1085.391, 1085.395, 1085.398, 1085.402, 1085.406, 1085.41, 1085.414, 1085.418, 1085.422, 1085.426], [4.78844635277511 3.1333248331567223 … -3.087774108506375 0.0032154107842656916; -0.47657443042814013 1.217125584744122 … 2.7853853316458257 0.0029729559817918234; … ; 6.441793000605642 5.047332343535441 … -0.6385598239986807 0.0008236881401195717; -27.729377847865074 -33.98805752002153 … -8.80813398237423 0.0018241535887852933;;;], Dict{Any, Any}(), 0×5 DataFrame
 Row │ id      start    length   description  channel 
     │ String  Float64  Float64  String       Int64   
─────┴────────────────────────────────────────────────, 23×9 DataFrame
 Row │ label   loc_radius  loc_theta  loc_x    loc_y    loc_z    loc_radius_sp ⋯
     │ String  Float64     Float64    Float64  Float64  Float64  Float64       ⋯
─────┼──────────────────────────────────────────────────────────────────────────
   1 │ Fp1           1.0       108.0    -0.31     0.95    -0.03            1.0 ⋯
   2 │ Fp2           1.0        72.0     0.31     0.95    -0.03            1.0
   3 │ F7            1.0       144.0    -0.81     0.59    -0.03            1.0
   4 │ F3            0.65      129.0    -0.55     0.67     0.5             1.0
   5 │ Fz            0.51       90.0     0.0      0.72     0.7             1.0 ⋯
   6 │ F4            0.65       51.0     0.55     0.67     0.5             1.0
   7 │ F8            1.0        36.0     0.81     0.59    -0.03            1.0
   8 │ T3            1.0       180.0    -1.0      0.0     -0.03            1.0
   9 │ C3            0.51      180.0    -0.72     0.0      0.7             1.0 ⋯
  10 │ Cz            0.0         0.0     0.0      0.0      1.0             1.0
  11 │ C4            0.51        0.0     0.72     0.0      0.7             1.0
  ⋮  │   ⋮         ⋮           ⋮         ⋮        ⋮        ⋮           ⋮       ⋱
  14 │ P3            0.65      231.0    -0.55    -0.67     0.5             1.0
  15 │ Pz            0.51      270.0     0.0     -0.72     0.7             1.0 ⋯
  16 │ P4            0.65      309.0     0.55    -0.67     0.5             1.0
  17 │ T6            1.0       324.0     0.81    -0.59    -0.03            1.0
  18 │ O1            1.0       252.0    -0.31    -0.95    -0.03            1.0
  19 │ O2            1.0       288.0     0.31    -0.95    -0.03            1.0 ⋯
  20 │ A1            1.0       192.0    -0.92    -0.23    -0.55            1.1
  21 │ A2            1.0       -12.0     0.92    -0.23    -0.55            1.1
  22 │ EOG1          1.01      149.0    -0.87     0.51    -0.37            1.0
  23 │ EOG2          1.01       31.0     0.87     0.51    -0.37            1.0 ⋯
                                                    3 columns and 2 rows omitted, ["reset_components(OBJ)", "filter(OBJ, ch=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24], fprototype=iirnotch, ftype=nothing, cutoff=50, order=8, rp=-1, rs=-1, dir=twopass, w=nothing)", "reset_components(OBJ)", "filter(OBJ, ch=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24], fprototype=fir, ftype=hp, cutoff=0.1, order=8, rp=-1, rs=-1, dir=twopass, w=nothing)", "reset_components(OBJ)", "trim(OBJ, seg=(1, 257), remove_epochs=false)", "reset_components(OBJ)", "trim(OBJ, seg=(1793, 2561), remove_epochs=false)", "reset_components(OBJ)", "trim(OBJ, seg=(1537, 2049), remove_epochs=false)"  …  "reset_components(OBJ)", "trim(OBJ, seg=(19137, 21313), remove_epochs=false)", "reset_components(OBJ)", "trim(OBJ, seg=(17409, 19585), remove_epochs=false)", "reset_components(OBJ)", "trim(OBJ, seg=(49537, 51265), remove_epochs=false)", "reset_components(OBJ)", "trim(OBJ, seg=(1, 5121), remove_epochs=false)", "reset_components(OBJ)", "ica_remove(OBJ, ch=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 19, 20, 21], ic_idx=1)"])

(!) ica_remove() reconstructs the source signal using only the specified component(s) and subtract the reconstructed signal from the source.

ICA: plotting

Plot the components:

p1 = NeuroAnalyzer.plot(eeg, ch="eeg", mono=true, seg=(0, 20), title="EEG channels")
p2 = NeuroAnalyzer.plot(eeg, ic, mono=true, seg=(0, 20), title="ICA components")
Plots.plot(p1, p2, layout=(2, 1))

Plot the components power spectrum:

p1 = plot_psd(eeg, ch="Fp1", title="original signal", mono=true)
p2 = plot_psd(eeg, ic, c_idx=2, title="ICA #02", mono=true)
eeg_new = ica_reconstruct(eeg, ic, ic_mw, ch="eeg", ic_idx=2, keep=true)
p3 = plot_psd(eeg_new, ch="Fp1", title="signal reconstructed from ICA #02", mono=true)
eeg_new = ica_reconstruct(eeg, ic, ic_mw, ch="eeg", ic_idx=2)
p4 = plot_psd(eeg_new, ch="Fp1", title="signal after removing ICA #02", mono=true)
Plots.plot(p1, p2, p3, p4, layout=(2, 2))

Plot the component spectrogram:

p1 = plot_spectrogram(eeg, seg=(0, 10), ch="Fp1", title="original signal")
p2 = plot_spectrogram(eeg, ic, c_idx=2, title="ICA #02")
eeg_new = ica_reconstruct(eeg, ic, ic_mw, ch="eeg", ic_idx=2, keep=true);
p3 = plot_spectrogram(eeg_new, ch="Fp1", title="signal reconstructed from ICA #02")
eeg_new = ica_reconstruct(eeg, ic, ic_mw, ch="eeg", ic_idx=2);
p4 = plot_spectrogram(eeg_new, ch="Fp1", title="signal after removing ICA #02")
Plots.plot(p1, p2, p3, p4, layout=(2, 2))