Convolution

Initialize NeuroAnalyzer 
using DSP
using NeuroAnalyzer
eeg = load("files/eeg.hdf")
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}(:transducers => ["", "", "", "", "", "", "", "", "", ""  …  "", "", "", "", "", "", "", "", "", ""], :epoch_id => "", :channel_order => [1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  15, 16, 19, 20, 21, 17, 18, 22, 23, 24], :channel_type => ["eeg", "eeg", "eeg", "eeg", "eeg", "eeg", "eeg", "eeg", "eeg", "eeg"  …  "eeg", "eeg", "ref", "ref", "eeg", "eeg", "eeg", "eog", "eog", "ecg"], :unit => ["μV", "μV", "μV", "μV", "μV", "μV", "μV", "μV", "μV", "μV"  …  "μV", "μV", "μV", "μV", "μV", "μV", "μV", "μV", "μV", "mV"], :file_size_mb => 52044.49, :recording_time => "12:32:53", :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"], :line_frequency => 50…), Dict(:name => "", :design => "", :notes => "")), ["filter_apply(OBJ, ch=[1], dir=twopass)", "filter_apply(OBJ, ch=[1], dir=twopass)", "filter_apply(OBJ, ch=[1], dir=twopass)", "filter_apply(OBJ, ch=[2], dir=twopass)", "filter_apply(OBJ, ch=[2], dir=twopass)", "filter_apply(OBJ, ch=[2], dir=twopass)", "filter_apply(OBJ, ch=[3], dir=twopass)", "filter_apply(OBJ, ch=[3], dir=twopass)", "filter_apply(OBJ, ch=[3], dir=twopass)", "filter_apply(OBJ, ch=[4], dir=twopass)"  …  "trim(OBJ, seg=(0, 1060), keep=true", "add_marker(OBJ, id=NA, start=62.0, len=1.0, value=DELETED, ch=0)", "trim(OBJ, seg=(62.0, 67.0), keep=false", "add_marker(OBJ, id=NA, start=74.0, len=1.0, value=DELETED, ch=0)", "trim(OBJ, seg=(74.0, 76.0), keep=false", "add_marker(OBJ, id=NA, start=119.0, len=1.0, value=DELETED, ch=0)", "trim(OBJ, seg=(119.0, 122.0), keep=false", "add_marker(OBJ, id=NA, start=348.0, len=1.0, value=DELETED, ch=0)", "trim(OBJ, seg=(348.0, 352.0), keep=false", "trim(OBJ, seg=(0, 1000), keep=true"], 0×5 DataFrame
 Row │ id      start    length   value   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 │ F3            0.65      129.0    -0.55     0.67     0.5             1.0
   4 │ F4            0.65       51.0     0.55     0.67     0.5             1.0
   5 │ C3            0.51      180.0    -0.72     0.0      0.7             1.0 ⋯
   6 │ C4            0.51        0.0     0.72     0.0      0.7             1.0
   7 │ P3            0.65      231.0    -0.55    -0.67     0.5             1.0
   8 │ P4            0.65      309.0     0.55    -0.67     0.5             1.0
   9 │ O1            1.0       252.0    -0.31    -0.95    -0.03            1.0 ⋯
  10 │ O2            1.0       288.0     0.31    -0.95    -0.03            1.0
  11 │ F7            1.0       144.0    -0.81     0.59    -0.03            1.0
  ⋮  │   ⋮         ⋮           ⋮         ⋮        ⋮        ⋮           ⋮       ⋱
  14 │ T4            1.0         0.0     1.0      0.0     -0.03            1.0
  15 │ T5            1.0       216.0    -0.81    -0.59    -0.03            1.0 ⋯
  16 │ T6            1.0       324.0     0.81    -0.59    -0.03            1.0
  17 │ A1            1.0       192.0    -0.92    -0.23    -0.55            1.1
  18 │ A2            1.0       -12.0     0.92    -0.23    -0.55            1.1
  19 │ Fz            0.51       90.0     0.0      0.72     0.7             1.0 ⋯
  20 │ Cz            0.0         0.0     0.0      0.0      1.0             1.0
  21 │ Pz            0.51      270.0     0.0     -0.72     0.7             1.0
  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, [0.0, 0.0039, 0.0078, 0.0117, 0.0156, 0.0195, 0.0234, 0.0273, 0.0312, 0.0352  …  999.9648, 999.9688, 999.9727, 999.9766, 999.9805, 999.9844, 999.9883, 999.9922, 999.9961, 1000.0], [0.0, 0.0039, 0.0078, 0.0117, 0.0156, 0.0195, 0.0234, 0.0273, 0.0312, 0.0352  …  999.9648, 999.9688, 999.9727, 999.9766, 999.9805, 999.9844, 999.9883, 999.9922, 999.9961, 1000.0], [-3.6187574610700066 -2.62108236965956 … 0.7186384231016238 -2.64668517785986; 0.45147186009170137 -0.614710955358406 … 1.926998524477149 -0.28398943425705087; … ; -0.9634647705466765 -2.076892851248747 … -14.603199667767157 -11.258329210808718; -5.436957946950784 -8.15171447987362 … -35.84200231407022 -37.465335411926226;;;])

Convolution is a mathematical operation that combines two functions to produce a third function. In signal processing, it involves sliding a kernel (or filter) along a signal and computing the dot product at each position.

This operation reveals how much the signal and kernel have in common at each point in time.

Mathematical Formulation

The convolution of a signal \(s\) and a kernel \(k\) is defined as:

\[ (s * k)[n] = \sum_{i=-\infty}^{\infty} s[i] \cdot k[n - i] \]

where:
\(s\) is the signal (e.g., EEG data),
\(k\) is the kernel (e.g., wavelet or filter),
\(n\) is the time index.

Sliding Dot Product

Convolution can be visualized as a sliding dot product between the signal and the kernel. The kernel is flipped (rotated) before sliding:

For a signal \(x = [3, 1, 4, 1, 5, 9]\) and a kernel \(w = [5, 7, 7]\), the convolution is computed as:

x = [3, 1, 4, 1, 5, 9]
w = [5, 7, 7]
conv(x, w)
8-element Vector{Int64}:
 15
 26
 48
 40
 60
 87
 98
 63

Kernel is rotated, so for w = [5, 7, 7] convolution is calculated with w = [7, 7, 5].

        [3   1   4   1   5   9]      n
[7   7   5]                          m

    [15                          ]   n + m - 1

        [3   1   4   1   5   9]
    [7   7   5]

    [15  26                      ]

        [3   1   4   1   5   9]
        [7   7   5]

    [15  26  48                  ]

The resulting convolved signal is longer than the original:

\[ l = N + M - 1 \]

where:
\(N\) is the signal length,
\(M\) is the kernel length.

Convolution vs Cross-Correlation vs Auto-Correlation

Group Delay

The group delay is the number of samples at the beginning of the signal before the convolution takes full effect:

\[ \tau_g = \frac{M - 1}{2} \]

where \(M\) is the length of the kernel (in samples).

Convolution Theorem

The convolution theorem states that convolution in the time domain is equivalent to multiplication in the frequency domain:

\[ \mathcal{F}\{s * k\} = \mathcal{F}\{s\} \cdot \mathcal{F}\{k\} \]

where \(\mathcal{F}\) is the Fourier transform.

Practical Use

  • In the time domain, convolution computes repeated dot products between the kernel and the signal.
  • In the frequency domain, convolution is equivalent to multiplying the Fourier transforms of the signal and kernel.

Circular Convolution

Circular convolution is a special case where the signal and kernel are treated as periodic. It is generally avoided in neuroscience due to artifacts, but it can be useful in specific applications like fast Fourier transform (FFT)-based convolution.

This tail is the copy of the initial segment of the convoluted signal.

This leads to time-domain aliasing. The solution is zero-padding to the next power of 2. For very long signals it is better to use block convolution.

Block convolution: when convoluting the signal in block, the tail of the convoluted segment must be added to the beginning of the next convoluted segment.

Both signal and kernel must be padded to the length: len(signal) + len(kernel) - 1.

After convolution padding has to be trimmed:

TRIM = 1/2 LENGTH(kernel) : SIGNAL : 1/2 LENGTH(kernel)

Kernel of odd length: convenient, but not required; easier trimming.

Role for EEG analysis:

  • isolating frequency-specific activity
  • localizing frequency-specific activity in time

Time-frequency analysis using wavelets:

  • create family of wavelets with different frequencies
  • perform convolution with each wavelet

Result of convolution shows when and to what extent the EEG data contains features that look like the wavelets. By using wavelets of different frequencies we can obtain time-frequency representation of the signal.

Time-Domain Convolution

mw = generate_morlet(32, 1, 32; complex = true)
tconv(eeg; ch = "all", kernel = mw)
24×256001×1 Array{ComplexF64, 3}:
[:, :, 1] =
  -11.5737+18.9615im    -14.9894+16.6061im   …   32.3805+15.6476im
   2.17636-8.87451im     4.02618-8.93077im      -23.9884+1.51382im
  -16.8092+25.859im     -21.8236+22.4243im      -5.16807+49.919im
  0.504024-21.4015im     4.37941-21.3793im       -16.672-10.4976im
  -8.13505+22.6609im    -13.7273+20.7157im      -32.1446+46.8684im
  -7.89086+1.07369im    -9.07813-0.679815im  …  -40.0392-18.2103im
   2.77519+34.0164im    -4.77584+34.0364im      -58.1395+37.6143im
   44.8043-27.1986im      48.789-19.1297im      -98.8087-68.297im
   11.7102+32.1714im     4.42741+34.0557im      -64.3975+45.4395im
   59.3269-25.3079im     63.2377-14.3781im      -126.152-59.3947im
   -22.344+30.8931im    -28.3483+26.4352im   …  -10.3867+36.1673im
  -9.49675-10.6891im    -8.51172-12.7524im      -33.8303-28.148im
   2.27057+38.5668im    -5.77996+38.9858im      -24.6318+60.9089im
   22.0033-15.083im      23.9812-11.4008im       -93.613-63.7162im
   10.3714+44.0341im     1.14767+46.0344im      -39.6121+58.4867im
   47.4179-19.9435im      50.529-11.3681im   …  -111.858-49.6278im
  -33.0614+28.9984im    -37.9925+21.4254im       26.5113+58.438im
   13.1864-56.9742im     24.7504-56.2326im      -94.7381-38.3085im
   2.93428+0.245609im    2.96876+0.19459im      -19.2529-5.36865im
  -18.3858+14.1662im    -22.3127+10.3257im       -55.153-3.93973im
   12.5142+16.1599im     7.87634+18.1267im   …   -59.339+20.6511im
   3.65736+26.1416im   -0.834819+26.35im         16.4062+59.8696im
 -0.991073-46.719im      8.49665-47.629im       -20.5661-45.5831im
   16.7388-35.5264im     24.0462-34.6527im      -118.286-314.691im

Frequency-Domain Convolution

eeg_new = deepcopy(eeg)
mw = generate_morlet(256, 1, 32, complex = true)
s_conv = fconv(eeg,
               ch = "all",
               kernel = mw)
eeg_new.data = abs.(s_conv)
24×256001×1 Array{Float64, 3}:
[:, :, 1] =
 0.148484  0.14873    0.148974   …  0.255401   0.255227   0.255048
 0.251746  0.252132   0.252512      0.180672   0.180325   0.179976
 0.348265  0.349352   0.350436      0.0692834  0.0692828  0.0692792
 0.06008   0.0602718  0.0604653     0.121341   0.121007   0.120672
 0.303222  0.304531   0.30584       0.193988   0.193315   0.192641
 0.201961  0.202985   0.204009   …  0.100249   0.099868   0.0994866
 0.434959  0.436257   0.437551      0.31272    0.311348   0.309977
 0.606436  0.608515   0.610592      0.109105   0.108102   0.107102
 0.417402  0.418902   0.4204        0.473529   0.471586   0.469643
 0.462833  0.464549   0.466262      0.356258   0.354836   0.353415
 0.288758  0.289812   0.290866   …  0.225283   0.224435   0.223586
 0.263537  0.264444   0.265348      0.111682   0.11125    0.110818
 0.398666  0.400231   0.401795      0.0864214  0.0865913  0.0867568
 0.400122  0.401616   0.40311       0.16367    0.163283   0.162894
 0.347284  0.348746   0.350209      0.211673   0.210429   0.209188
 0.39139   0.392994   0.394598   …  0.335089   0.334489   0.333883
 0.757163  0.759977   0.762789      0.346765   0.345492   0.344218
 0.564593  0.567235   0.56988       0.429435   0.42755    0.425665
 0.167516  0.168152   0.168787      0.0681442  0.0679461  0.0677473
 0.288675  0.290004   0.291335      0.0846258  0.0841902  0.0837557
 0.524455  0.526185   0.527911   …  0.236233   0.235183   0.234134
 0.212829  0.21291    0.212983      0.418041   0.417098   0.416149
 0.111027  0.112028   0.113037      0.120073   0.120634   0.12119
 2.24128   2.24905    2.2568        2.37307    2.36412    2.35517

Tip: Since kernel is complex-valued, the returned signal is also complex-valued. To replace the original signal, its absolute values must be used:

\[ \vert z \vert = \sqrt{\Re^2 + \Im^2} \]

Denoising using Wiener Deconvolution

denoise_wien(eeg,
             ch = "all")
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 => "")), ["filter_apply(OBJ, ch=[1], dir=twopass)", "filter_apply(OBJ, ch=[1], dir=twopass)", "filter_apply(OBJ, ch=[1], dir=twopass)", "filter_apply(OBJ, ch=[2], dir=twopass)", "filter_apply(OBJ, ch=[2], dir=twopass)", "filter_apply(OBJ, ch=[2], dir=twopass)", "filter_apply(OBJ, ch=[3], dir=twopass)", "filter_apply(OBJ, ch=[3], dir=twopass)", "filter_apply(OBJ, ch=[3], dir=twopass)", "filter_apply(OBJ, ch=[4], dir=twopass)"  …  "add_marker(OBJ, id=NA, start=62.0, len=1.0, value=DELETED, ch=0)", "trim(OBJ, seg=(62.0, 67.0), keep=false", "add_marker(OBJ, id=NA, start=74.0, len=1.0, value=DELETED, ch=0)", "trim(OBJ, seg=(74.0, 76.0), keep=false", "add_marker(OBJ, id=NA, start=119.0, len=1.0, value=DELETED, ch=0)", "trim(OBJ, seg=(119.0, 122.0), keep=false", "add_marker(OBJ, id=NA, start=348.0, len=1.0, value=DELETED, ch=0)", "trim(OBJ, seg=(348.0, 352.0), keep=false", "trim(OBJ, seg=(0, 1000), keep=true", "denoise_wien(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])"], 0×5 DataFrame
 Row │ id      start    length   value   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 │ F3            0.65      129.0    -0.55     0.67     0.5             1.0
   4 │ F4            0.65       51.0     0.55     0.67     0.5             1.0
   5 │ C3            0.51      180.0    -0.72     0.0      0.7             1.0 ⋯
   6 │ C4            0.51        0.0     0.72     0.0      0.7             1.0
   7 │ P3            0.65      231.0    -0.55    -0.67     0.5             1.0
   8 │ P4            0.65      309.0     0.55    -0.67     0.5             1.0
   9 │ O1            1.0       252.0    -0.31    -0.95    -0.03            1.0 ⋯
  10 │ O2            1.0       288.0     0.31    -0.95    -0.03            1.0
  11 │ F7            1.0       144.0    -0.81     0.59    -0.03            1.0
  ⋮  │   ⋮         ⋮           ⋮         ⋮        ⋮        ⋮           ⋮       ⋱
  14 │ T4            1.0         0.0     1.0      0.0     -0.03            1.0
  15 │ T5            1.0       216.0    -0.81    -0.59    -0.03            1.0 ⋯
  16 │ T6            1.0       324.0     0.81    -0.59    -0.03            1.0
  17 │ A1            1.0       192.0    -0.92    -0.23    -0.55            1.1
  18 │ A2            1.0       -12.0     0.92    -0.23    -0.55            1.1
  19 │ Fz            0.51       90.0     0.0      0.72     0.7             1.0 ⋯
  20 │ Cz            0.0         0.0     0.0      0.0      1.0             1.0
  21 │ Pz            0.51      270.0     0.0     -0.72     0.7             1.0
  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, [0.0, 0.0039, 0.0078, 0.0117, 0.0156, 0.0195, 0.0234, 0.0273, 0.0312, 0.0352  …  999.9648, 999.9688, 999.9727, 999.9766, 999.9805, 999.9844, 999.9883, 999.9922, 999.9961, 1000.0], [0.0, 0.0039, 0.0078, 0.0117, 0.0156, 0.0195, 0.0234, 0.0273, 0.0312, 0.0352  …  999.9648, 999.9688, 999.9727, 999.9766, 999.9805, 999.9844, 999.9883, 999.9922, 999.9961, 1000.0], [-3.6140631348209347 -2.649020174736541 … 0.7256278458920968 -2.657624532869198; 0.4176420853081487 -0.6038876857674688 … 1.9054175068021217 -0.3031541488360958; … ; -1.0121668004674111 -2.0381610531493135 … -14.602564831151941 -11.24571134072962; -5.428324222122386 -8.081085025292287 … -35.809125140409066 -37.45546465210564;;;])