Find EMG onset time using EMG & smoothed EMG signal
[zx] = vb_get_emg_onset_time(y ,yave, y0, z0, parm)
[zx, ixx ,jxx, kx] = vb_get_emg_onset_time(y ,yave, y0, z0, parm)
--- Input
y : EMG signal
yave : smoothed EMG by moving average
y0 : threshold for abs(y)
z0 : smoothed EMG threshold value
parm : parameter structure
parm.fsamp : sampling frequency [Hz]
--- Optional parameter for EMG onset
parm.t_event : minimum distance from previous onset event [150 ms]
--- Usually following parameters need not be changed
parm.t_period : minimum period that Smoothed EMG > threshold [50 ms]
--- Condition for EMG onset
1. Smoothed EMG (yave) should exceed 'z0' more than 't_period' length
2. Distance from previous onset should be larger than 't_event'
--- Output
zx : extracted onset time index
ixx : abs(EMG) onset time index
jxx : smoothed EMG onset
kx : start & end of slope estimation region at jx
2009-08-05 Masa-aki Sato
2011-12-18 Masa-aki Sato
Conditions and parameters for onset are changed
2012-2-18 Masa-aki Sato
Use TKE operator EMG filter
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function [zx ,ix_y, ix_z, kx] = vb_get_emg_onset_time(y,yave,y0,z0,parm) 0002 % Find EMG onset time using EMG & smoothed EMG signal 0003 % [zx] = vb_get_emg_onset_time(y ,yave, y0, z0, parm) 0004 % [zx, ixx ,jxx, kx] = vb_get_emg_onset_time(y ,yave, y0, z0, parm) 0005 % --- Input 0006 % y : EMG signal 0007 % yave : smoothed EMG by moving average 0008 % y0 : threshold for abs(y) 0009 % z0 : smoothed EMG threshold value 0010 % 0011 % parm : parameter structure 0012 % parm.fsamp : sampling frequency [Hz] 0013 % --- Optional parameter for EMG onset 0014 % parm.t_event : minimum distance from previous onset event [150 ms] 0015 % --- Usually following parameters need not be changed 0016 % parm.t_period : minimum period that Smoothed EMG > threshold [50 ms] 0017 % --- Condition for EMG onset 0018 % 1. Smoothed EMG (yave) should exceed 'z0' more than 't_period' length 0019 % 2. Distance from previous onset should be larger than 't_event' 0020 % 0021 % --- Output 0022 % zx : extracted onset time index 0023 % ixx : abs(EMG) onset time index 0024 % jxx : smoothed EMG onset 0025 % kx : start & end of slope estimation region at jx 0026 % 0027 % 2009-08-05 Masa-aki Sato 0028 % 2011-12-18 Masa-aki Sato 0029 % Conditions and parameters for onset are changed 0030 % 2012-2-18 Masa-aki Sato 0031 % Use TKE operator EMG filter 0032 % 0033 % Copyright (C) 2011, ATR All Rights Reserved. 0034 % License : New BSD License(see VBMEG_LICENSE.txt) 0035 0036 % Default value 0037 t_slope = 25; % 25 ms 0038 t_period = 50; % 50 ms 0039 t_event = 500; % 500 ms 0040 0041 if isfield(parm,'t_slope'), t_slope = parm.t_slope ;end 0042 if isfield(parm,'t_period'), t_period = parm.t_period ;end 0043 if isfield(parm,'t_event'), t_event = parm.t_event ;end 0044 0045 % sample number in [msec] 0046 ds = parm.fsamp / 1000; 0047 0048 % msec -> sample number 0049 t_slope = round(t_slope * ds); 0050 t_period = round(t_period * ds); 0051 t_event = round(t_event * ds); 0052 0053 % EMG threshold 0054 y = abs(y(:)'); 0055 yave = yave(:)'; 0056 T = length(y); 0057 0058 % Extract smoothed EMG onset (time at which yave exceed z0) 0059 % 'yave' value shoud be larger than 'z0' for 't_period' 0060 ix = vb_trigger_onset(yave,z0,t_period,t_period); 0061 0062 % time length from previous onset should be larger than 't_event' 0063 [ixx, jj] = vb_onset_check(ix,t_event); 0064 0065 % output for matched 'ix' and 'jx' 0066 ix_y = {}; 0067 ix_z = {ix; ixx}; 0068 0069 % onset candidate 0070 zx = ixx; 0071 kx = []; 0072 0073 %return 0074 % --- No zero crossing Estimation --- 0075 0076 if isempty(t_slope), return; end; 0077 if t_slope == 0, return; end; 0078 0079 % --- Estimate zero crossing --- 0080 % by linear fit to smoothed EMG near onset 0081 0082 % number of matched onset 0083 NN = length(ixx); 0084 0085 % --- Positive (non-negative) slope region for yave 0086 T = length(yave); 0087 0088 ix_sp = find( [0 diff(yave)] >= 0 ); 0089 slope = zeros(1, T ); 0090 slope(ix_sp) = 1; 0091 0092 % [0 slope] : [ 0 0 1 1 1 0 0 1 1 1 1 0 0 0 1 1 1] 0093 % diff : [ 0 1 0 0 -1 0 1 0 0 0 -1 0 0 1 0 0] 0094 0095 % start of positive slope 0096 k1 = find( diff([0 slope]) > 0 ); 0097 % end of positive slope 0098 k2 = find( diff([0 slope]) < 0 ); 0099 0100 N1 = min(length(k1),length(k2)); 0101 0102 % start & end of positive slope region 0103 k1 = k1(1:N1); 0104 k2 = k2(1:N1); 0105 0106 % start & end of slope estimation region 0107 kx = zeros(2,NN); 0108 0109 % Estimate zero crossing 0110 % by linear fit to smoothed EMG near onset 0111 for n=1:NN 0112 % onset candidate for smoothed EMG 0113 j0 = ixx(n); 0114 0115 % find positive slope region that include 'j0' 0116 kk = find( k1 < j0 & k2 > j0 ); 0117 0118 if isempty(kk), 0119 j1 = j0 - t_slope; 0120 j2 = j0 + t_slope; 0121 else 0122 j1 = min( k1(kk(1)) , j0 - t_slope); 0123 j2 = max( k2(kk(1)) , j0 + t_slope); 0124 end 0125 0126 j1 = max(j1,1); 0127 j2 = min(j2,T); 0128 0129 % linear fit to positive slope region 0130 % z - <z> = a * (x - <x>) 0131 % z = a * (x - <x>) + <z> = 0 0132 % x = <x> - <z>/a 0133 x = j1:j2; 0134 z = yave(x); 0135 xm = mean(x); 0136 zm = mean(z); 0137 a = sum((z-zm).*(x-xm))/sum((x-xm).^2); 0138 zx(n) = xm - zm/a; 0139 kx(:,n) = [j1 ; j2]; 0140 end 0141 0142 % check bad condition for linear fit 0143 ii = find( abs(ixx - zx) > t_slope ); 0144 0145 zx(ii) = ixx(ii); 0146 0147 return