-
Notifications
You must be signed in to change notification settings - Fork 6
/
run_decomposition.m
119 lines (115 loc) · 5.64 KB
/
run_decomposition.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
function motor_unit = run_decomposition(varargin)
%RUN_DECOMPOSITION motor-unit decomposition on hdEMG datasets
%
% This function and the helper files are mainly a reorganization of the
% code accompanied with Hyser Dataset by Jian et. al. <a href=
% 'https://doi.org/10.1109/TNSRE.2021.3082551'>Hyser Dataset</a>
% The original code and the dataset is also <a href=
% 'https://doi.org/10.13026/ym7v-bh53'>available at PhysioNet</a>
%
% Requirements:
% Matlab R2021b+, Statistics and Machine Learning Toolbox, DSP
% Toolbox, and Communications Toolbox (to inject white noise).
% This function uses tables, with features available from
% R2017b+.
%
% INPUT:
% Inputs are name pairs:
%
% 'data' : (The path to the) hdEMG data. if 'data' is pointing to the
% lcation of the data file, the data file must be a MAT array.
% Default is the sample file included in the toolbox.
%
% 'data_mode' : EMG can be recorded in the 'monopolar' or 'bipolar'
% mode. Default = 'monopolar'
%
% 'frq','sampling_frequency' : Sampling frequency of the data, default = 2048 Hz
%
% 'arr', 'array_shape' : shape of the arry, as a 1 x 2 vevtor. Default =
% [8, 8]. Note, setting 'arr', [13, 5] will throw an error; please use [8, 8] instead.
% This is because the computing functions are not spatially sensitve. Later on,
% and for array plots, the array shape will determine the plot shape
%
% 'R', 'extension_parameter' : The number of times to repeat the data
% blocks, see the Hyser paper for more detail. Default = 4
%
% 'M', 'max_iter' : Maximum iterations for the (FAST) ICA decompsition.
% Default = 300
%
% 'whiten_flag' : Whether to whiten the data prior to the ICA. Default = 1
%
% 'SNR', 'inject_noise' : Adding white noise to the EMG mixutre. Uses
% Communication Toolbox AWGN fucntion. Default = Inf for not injecting
% any artificial noise.
%
% 'SIL_thresh' : The silhouette threshold to detect the good motor units.
% Default = 0.6
%
% 'save_path' : The path that the files should be saved there. The
% function does not create the path, rather uses it. Default is the
% is the 'sample' path of the toolbox.
%
% 'save_flag' : Whether the files are saved or not, default is 0, so it is
% NOT saving your output.
%
% OUTPUT:
% 'motor_unit': The structure inclduing the follwing fields:
% spike_train
% waveform
% ica_weight
% whiten_matrix
% SIL
%
% EXAMPLE:
% sample_motorUnit = RUN_DECOMPOSITION()
%
% REV:
% v0 @ 09/12/2022
%
%
% Copyright (c) 2022 Seyed Yahya Shirazi, [email protected]
%% initialize
fs = filesep; % a shortcut to the filesep
addpath('./helper_functions')
opts = arg_define(varargin, ...
arg({'data'},['sample_data' fs 'sample1'] ,[],'The file is also from the Hyser dataset'), ...
arg({'data_mode'}, 'monopolar',['monopolar','bipolar'],'Data mode of the hd-EMG data'), ...
arg({'frq','sampling_frequency'}, 2048,[],'Sampling frequency of the imported dataset'), ...
arg({'array_shape','arr'}, [8, 8],[],'Array shape as a 2D vector'), ...
arg({'R', 'extension_parameter'}, 4,[],'The number of times to repeat the data blocks'), ...
arg({'M', 'max_sources'}, 300,[],'Maximum number of sources being decomposed by (FAST) ICA.'), ...
arg({'whiten_flag'}, 1,[0,1],'Whether to whiten the data prior to the ICA.'), ...
arg({'SNR', 'inject_noise'}, Inf,[],'WAdding white noise to the EMG mixutre.'), ...
arg({'SIL_thresh'}, 0.6,[],'The silhouette threshold to detect the good motor units.'), ...
arg({'output_file'},['sample_data' fs 'sample1_decomposed'] ,[],'The path to the output file.'), ...
arg({'save_flag'}, 1,[0,1],'Whether the files are saved or not.'),...
arg({'plot_spikeTrain', 'plot_fig'}, 1,[0,1],'Whether to plot the resulting spike trains.'),...
arg({'load_ICA'}, 0,[0,1],'Whether to load precomputed ICA results for debugging.'));
data = opts.data;
if ischar(data)
emg_file = load(data); % Based on the strcutre of the exproted MAT files from OTBiolab+ v1.5.8
data = emg_file.Data{1}; % Assuming that the file ONLY contains one hdEMG array!
opts.frq = emg_file.SamplingFrequency;
end
max_iter = 200;
%% run the decomposition
[extended_emg,~] = emg_preprocess(data,'R',opts.R,'WhitenFlag',opts.whiten_flag,'SNR',opts.SNR, 'array_shape', opts.array_shape);
if ~opts.load_ICA
[uncleaned_source,B,uncleaned_spkieTrain,score] = run_ICA(extended_emg, opts.M, max_iter); % B is the unmixing matrix
else
warning("Loading ICA results from a saved file. Change the 'load_ICA' flag if you want to run ICA.")
ica_results = load([opts.data '_ica_results']);
uncleaned_source = ica_results.uncleaned_source; uncleaned_spkieTrain = ica_results.uncleaned_spkieTrain;
score = ica_results.score; B = ica_results.B;
end
[spike_train, source, good_idx] = remove_motorUnit_duplicates(uncleaned_spkieTrain, uncleaned_source, opts.frq);
silhouette_score = score(good_idx);
%% save the results
% motor_unit.raw.B = B; motor_unit.raw.extended_emg = extended_emg; % this creates a very large file thay would not fit on github.
% motor_unit.raw.uncleaned_source = uncleaned_source; motor_unit.raw.score = score;
motor_unit.spike_train = spike_train; motor_unit.source = source; motor_unit.good_idx = good_idx;
motor_unit.silhouette_score = silhouette_score;
if opts.save_flag, save(opts.output_file,"-struct","motor_unit"); end
%% plot the motor units in a nice way
minScore_toPlot = 0.9; % The minimum silhouette score of the motor units to be included in the plot
if opts.plot_spikeTrain, plot_spikeTrain(spike_train,opts.frq,silhouette_score,minScore_toPlot); end