vb_test_emg_threshold

 Script for check EMG threshold for vb_job_trial_onset
 --- Input
 parm.data_file : Data file name         [string]
 parm.trig_file : Trial onset file name  [string]
 parm.status_ch : status channel name    [string]
 parm.trig_type = 'integer' or 'analog' or 'voice' or 'emg'
 parm.slope = 'const_start' or 'const_end'    if trig_type = 'integer'
              'low_to_high' or 'high_to_low'  if trig_type = 'analog'
              No meaning                      if trig_type = 'voice','emg'
 parm.condition : string describing condition [string or cell array]
 parm.status_level : status level      [1 x Ncomdition]
 parm.Pretrigger_ms : Pretrigger period   [msec]
 parm.Posttrigger_ms: Posttrigger period  [msec]

0001 % Script for check EMG threshold for vb_job_trial_onset
0002 % --- Input
0003 % parm.data_file : Data file name         [string]
0004 % parm.trig_file : Trial onset file name  [string]
0005 % parm.status_ch : status channel name    [string]
0006 % parm.trig_type = 'integer' or 'analog' or 'voice' or 'emg'
0007 % parm.slope = 'const_start' or 'const_end'    if trig_type = 'integer'
0008 %              'low_to_high' or 'high_to_low'  if trig_type = 'analog'
0009 %              No meaning                      if trig_type = 'voice','emg'
0010 % parm.condition : string describing condition [string or cell array]
0011 % parm.status_level : status level      [1 x Ncomdition]
0012 % parm.Pretrigger_ms : Pretrigger period   [msec]
0013 % parm.Posttrigger_ms: Posttrigger period  [msec]
0014 
0015 % --- Optional parameter for EMG onset
0016 % parm.t_event  : minimum distance from previous onset event [150 ms]
0017 %                 distance from previous onset should be larger than t_event
0018 % parm.p_val : P-value corresponding to the threshold for [EMG, smooth(EMG)]
0019 %               [0.0001, 0.0005] or [0.0001, 0.001]
0020 %     cumulative histgram is used to determine threshold from P-value
0021 %
0022 % --- Usually following parameters need not be changed
0023 % parm.hist_mode : histgram mode [1]
0024 %                = 1: Estimate threshold by gamma distribution approximation
0025 %                = 0: Estimate threshold by histgram
0026 % parm.t_smooth : moving average window length               [25 ms]
0027 % parm.t_slope  : slope estimation period near onset         [25 ms]
0028 %                 if t_slope==0, zero cross point is not calculated
0029 % parm.t_peak   : peak evaluation period                     [100 ms]
0030 % peak_val : EMG value should exceed peak_val within 't_peak' after onset
0031 %          = mean(peak value > threshold) if hist_mode = 1
0032 %          = (max(y) * status_level)      if hist_mode = 0 or 2
0033 % --- Condition for EMG onset
0034 % 1. distance from previous onset should be larger than t_event
0035 % 2. distance between EMG & smoothed EMG onset should be smaller than t_event
0036 % 3. EMG value should exceed peak_val within t_peak after onset
0037 % 4. zero cross point is estimated
0038 %    by linear fitting around smoothed EMG threshold point
0039 %
0040 % 2009-09-05 Masa-aki Sato
0041 %
0042 % Copyright (C) 2011, ATR All Rights Reserved.
0043 % License : New BSD License(see VBMEG_LICENSE.txt)
0044 
0045 clear all
0046 
0047 data_mode = 1;
0048 plot_mode = 1;
0049     % 0: Plot EMG and onset + EMG histgram
0050     % 1: Plot EMG and onset and threshold
0051     % 2: Plot EMG and onset for each onset in subplot
0052     
0053     % 3: Plot onset interval
0054     % 4: Plot histgram
0055     % 5: Plot histgram detail
0056 tlim = [415 420];
0057 
0058 % --- Please change directry and file name
0059 % MEG-mat root directry
0060 switch    data_mode
0061 case    1
0062     root_dir  = [getenv('MATHOME') '/Finger_move/RH070222'];
0063     data_file = 'RH070222a.meg.mat';
0064     status_ch = '226'; % channel name for EMG signal：226,227
0065 
0066     loadspec = [];
0067     loadspec.ChannelName = {status_ch};
0068     status = vb_load_meg_data([root_dir '/' data_file], loadspec);
0069     
0070     info = vb_load_meg_info([root_dir '/' data_file]);
0071     parm.fsamp = info.SampleFreq;
0072     
0073     % setting for plot
0074     Ntrial = 10; % onset number in one subplot
0075     t_step = 200; % Xtick step [msec]
0076     Tpre   = 100; % pre time for onset plot [msec]
0077     Tpost  = 500; % post time for onset plot [msec]
0078     Tend   = 5000;
0079 case    2
0080     root_dir  = './test'; 
0081     data_file = 'emg.mat';
0082     load( [root_dir '/' data_file] ,'status');
0083     parm.fsamp = 1000;  % [Hz]
0084     Ntrial = 25; % onset number in one subplot
0085     t_step = 500;
0086     Tpre   = 500; 
0087     Tpost  = 800; 
0088     Tend   = 5000;
0089 end
0090 
0091 % parm.p_val : P-value for [EMG, smooth(EMG)] threshold
0092 parm.p_val = [0.0001 0.0005];
0093 %parm.p_val = [0.0001 0.001];
0094 
0095 % --- Set minimum distance from previous onset event [ms]
0096 switch    data_mode
0097 case    1
0098     parm.t_event  = 150;
0099 case    2
0100     parm.t_event  = 300;
0101 end
0102 
0103 % --- Usually following parameters need not be changed
0104 parm.hist_mode = 1;
0105 
0106 % parm.t_smooth : moving average window length        [25 ms]
0107 % parm.t_slope  : slope estimation period near onset  [25 ms]
0108 % parm.t_peak   : peak evaluation period              [100 ms]
0109 parm.t_smooth = 25 ;
0110 parm.t_slope  = 25 ;
0111 parm.t_peak   = 100;
0112 
0113 % EMG signal
0114 y  = status(:);
0115 
0116 dt = 1/parm.fsamp; % [sec]
0117 T  = length(y);
0118 t  = (1:T)*dt;
0119 
0120 % EMG onset
0121 % z    : smoothed EMG by moving average
0122 % zd   : smoothed EMG threshold value
0123 % yd   : threshold for abs(y)
0124 %
0125 % ix  : extracted onset time index
0126 % ix_y  : abs(EMG) onset time index
0127 % ix_z  : smoothed EMG onset
0128 % ix_s  : start & end of slope estimation region at ix_z
0129 
0130 [ix ,zd, z, ix_y, yd, ix_z, ix_s] = vb_get_emg_onset(y,parm);
0131 
0132 % onset time
0133 jz  = ix * dt;
0134 j_z = ix_z * dt;
0135 
0136 ymax  = max(abs(y));
0137 % X-tick
0138 xtick = 0:(t_step/1000):ceil(t(end));
0139 % onset display time window
0140 t1 = jz - Tpre/1000;
0141 t2 = jz + Tpost/1000;
0142 
0143 Nz  = length(ix);
0144 nlist = 1:Nz;
0145 %nlist = [10:15];
0146 
0147 switch    plot_mode
0148 case    0
0149     % Plot EMG and threshold
0150     
0151     % active check threshold
0152     yth = vb_get_peak_threshold(abs(y),yd);
0153     yd2 = (yd + mean(yth))/2;
0154     %yd2 = max(y) * parm.status_level;
0155 
0156     figure;
0157     subplot    211
0158     % Plot EMG
0159     plot(t,abs(y),'--')
0160     hold on
0161     % estimated onset
0162     plot([jz; jz], [zeros(1,Nz); repmat(1.1*ymax, 1,Nz)], 'r--');
0163     % onset threshold
0164     plot([0; t(end)], [yd; yd], 'r-')
0165     plot([0; t(end)], [yd2; yd2], 'b--')
0166     title('EMG threshold and onset')
0167     
0168     subplot    212
0169     % Plot smoothed EMG
0170     plot(t,z,'c-')
0171     hold on
0172 
0173     % estimated onset
0174     plot([jz; jz], [zeros(1,Nz); repmat(ymax, 1,Nz)], 'r-');
0175     % onset threshold
0176     plot([0; t(end)], [zd; zd], 'r-')
0177     title('Smoothed EMG threshold and onset')
0178     
0179     figure;
0180     
0181     % EMG histgram
0182     Nbin = 1000;
0183     [hy ,ylist] = hist(abs(y),Nbin);
0184     [hz ,zlist] = hist(z,Nbin);
0185 
0186     subplot 221
0187     % EMG histgram
0188     hist(abs(y),Nbin);
0189     hold on
0190     plot([yd; yd], [0; max(hy)], 'r-')
0191     plot([yd2; yd2], [0; max(hy)], 'b--')
0192     
0193     title('EMG histogram')
0194     
0195     subplot 222
0196     % Smoothed EMG histgram
0197     hist(z,Nbin);
0198     hold on
0199     plot([zd; zd], [0; max(hz)], 'r-')
0200     title('Smoothed EMG histogram')
0201     
0202     subplot 223
0203     % Onset interval
0204     jzd = diff(jz);
0205     plot(jzd,'o')
0206     title('Onset interval')
0207 case    1
0208     % Plot EMG and threshold
0209     % active check threshold
0210     yth = vb_get_peak_threshold(abs(y),yd);
0211     yd2 = (yd + mean(yth))/2;
0212     %yd2 = max(y) * parm.status_level;
0213     Nfig = ceil(Nz/Ntrial);
0214     t_wd = Tend/1000;
0215     
0216     for n=1:Nfig
0217         n1 = 1+Ntrial*(n-1);
0218         n2 = min(Nz,n1+Ntrial-1);
0219         figure;
0220         subplot    211
0221         % Plot EMG
0222         plot(t,abs(y),'--')
0223         hold on
0224         
0225         % estimated onset
0226         plot([jz; jz], [zeros(1,Nz); repmat(1.1*ymax, 1,Nz)], 'r--');
0227         % onset threshold
0228         plot([0; t(end)], [yd; yd], 'r-')
0229         plot([0; t(end)], [yd2; yd2], 'b--')
0230         title('EMG threshold and onset')
0231         
0232         xlim([jz(n1)-t_wd jz(n2)+t_wd])
0233         
0234         subplot    212
0235         % Plot smoothed EMG
0236         plot(t,z,'c-')
0237         hold on
0238     
0239         % estimated onset
0240         plot([jz; jz], [zeros(1,Nz); repmat(ymax, 1,Nz)], 'r--');
0241         % onset threshold
0242         plot([0; t(end)], [zd; zd], 'r-')
0243         title('Smoothed EMG threshold and onset')
0244         xlim([jz(n1)-t_wd jz(n2)+t_wd])
0245     end
0246 case    {2 , 20}
0247     % Plot EMG and onset in subplot
0248     nfig=0; npng = 0;
0249     NX = 3; NY = 4;
0250 
0251     for n=nlist
0252         nfig = nfig + 1;
0253         if nfig > NX*NY, 
0254             npng = npng + 1;
0255             nstr = num2str(npng);
0256             if plot_mode==2
0257                 print(gcf, '-dpng', ['temp_' nstr '.png']);
0258                 close
0259             end
0260             figure; nfig=1; 
0261         end
0262         subplot(NY,NX,nfig)
0263         
0264         plot(t,abs(y),':')
0265         hold on
0266         plot(t,z,'c-')
0267         
0268         % threshold
0269         plot([0; t(end)], [yd; yd], 'b-')
0270         plot([0; t(end)], [zd; zd], 'b--')
0271         
0272         % estimated onset
0273         plot([jz; jz], [zeros(1,Nz); repmat(ymax, 1,Nz)], 'r--');
0274         plot(j_z, z(ix_z), 'ro');
0275         
0276         % slope estimation period
0277         if ix_s(1,n) > 0
0278             %plot(t(ix_s(1,n):ix_s(2,n)) , z(ix_s(1,n):ix_s(2,n)), 'r-')
0279             plot(t(ix_s(1,n):ix_s(2,n)) , z(ix_s(1,n):ix_s(2,n)),...
0280                  'r-','LineWidth',2)
0281         end
0282         
0283         xlim([t1(n)  t2(n)])
0284         ylim([0 ymax])
0285         set(gca,'XTick',xtick);
0286         title(sprintf('n=%d',n))
0287         
0288         
0289     end
0290 case    3
0291     % Onset interval
0292     jzd = diff(jz);
0293     plot(jzd,'o')
0294     title('Onset interval')
0295     
0296 case    4
0297     Nbin = 1000;
0298     [hy ,ylist] = hist(abs(y),Nbin);
0299     [hz ,zlist] = hist(z,Nbin);
0300     yth = vb_get_peak_threshold(abs(y),yd);
0301     yd2 = (yd + mean(yth))/2;
0302     %yd2 = max(y) * parm.status_level;
0303 
0304     subplot 221
0305     hist(abs(y),Nbin);
0306     hold on
0307     plot([yd; yd], [0; max(hy)], 'r-')
0308     
0309     title('EMG histogram')
0310     
0311     subplot 222
0312     hist(z,Nbin);
0313     hold on
0314     plot([zd; zd], [0; max(hz)], 'r-')
0315     title('smoothed EMG histogram')
0316 
0317     subplot 223
0318     plot(t,abs(y))
0319     hold on
0320     plot([t(1) t(end)], [yd; yd],'r-')
0321     plot([t(1) t(end)], [yd2; yd2],'r--')
0322     title('EMG signal')
0323 
0324     subplot 224
0325     plot(t,z)
0326     hold on
0327     plot([t(1) t(end)], [zd; zd],'r-')
0328     title('smoothed EMG signal')
0329 case    5
0330     status_level = 1/4;
0331     Nbin = min(5000,fix(T/10));
0332     yd2 = max(y) * status_level;
0333     
0334     y = abs(y);
0335     [hy ,ylist] = hist(y,Nbin);
0336     [hz ,zlist] = hist(z,Nbin);
0337     
0338     [yth, A ,ym, yy] = vb_gamma_dist_param(y, parm.p_val(1));
0339     py = vb_gamma_dist_prob(ylist,hy,A);
0340     yd3 = vb_get_peak_threshold(y,yth)
0341 
0342     [zth, A ,zm, zz] = vb_gamma_dist_param(z, parm.p_val(2));
0343     pz = vb_gamma_dist_prob(zlist,hz,A);
0344 
0345     subplot 211
0346     hist(abs(y),Nbin);
0347     hold on
0348     plot([yd; yd], [0; max(hy)], 'r-')
0349     plot([yth; yth], [0; max(hy)], 'b-')
0350 
0351     plot([yd2; yd2], [0; max(hy)], 'r:')
0352     plot([yd3(1); yd3(1)], [0; max(hy)], 'b:')
0353     plot([yd3(2); yd3(2)], [0; max(hy)], 'b:')
0354     
0355     plot([ym; ym], [0; max(hy)], 'r--')
0356     plot([ym+yy; ym+yy], [0; max(hy)], 'r--')
0357     plot([ym-yy; ym-yy], [0; max(hy)], 'r--')
0358     plot(ylist,py,'-r')
0359     
0360     xlim([0 1.1*max([yd,yth,yd2,yd3])])
0361     
0362     title('EMG histogram')
0363     
0364     subplot 212
0365     hist(z,Nbin);
0366     hold on
0367     plot([zd; zd], [0; max(hz)], 'r-')
0368     plot([zth; zth], [0; max(hz)], 'b-')
0369     plot([zm; zm], [0; max(hz)], 'r--')
0370     plot([zm+zz; zm+zz], [0; max(hz)], 'r--')
0371     plot([zm-zz; zm-zz], [0; max(hz)], 'r--')
0372     plot(zlist,pz,'-r')
0373     
0374     xlim([0 1.1*max(zd,zth)])
0375     title('smoothed EMG histogram')
0376 case    10
0377     % Plot EMG and threshold
0378     % active check threshold
0379     yth = vb_get_peak_threshold(abs(y),yd);
0380     yd2 = (yd + mean(yth))/2;
0381     figure;
0382     subplot    211
0383     % Plot EMG
0384     plot(t,abs(y),'--')
0385     hold on
0386     
0387     % estimated onset
0388     plot([jz; jz], [zeros(1,Nz); repmat(1.1*ymax, 1,Nz)], 'r--');
0389     % onset threshold
0390     plot([0; t(end)], [yd; yd], 'r-')
0391     plot([0; t(end)], [yd2; yd2], 'b--')
0392     title('EMG threshold and onset')
0393     
0394     xlim(tlim)
0395     
0396     subplot    212
0397     % Plot smoothed EMG
0398     plot(t,z,'c-')
0399     hold on
0400     
0401     % estimated onset
0402     plot([jz; jz], [zeros(1,Nz); repmat(ymax, 1,Nz)], 'r--');
0403     % onset threshold
0404     plot([0; t(end)], [zd; zd], 'r-')
0405     title('Smoothed EMG threshold and onset')
0406     xlim(tlim)
0407 end