vb_get_emg_threshold

 Estimate threshold for smoothed EMG power
    [y0, z0, z] = vb_get_emg_threshold(y,parm)
    
 --- Input
 y             : EMG signal
 
 parm.fsamp    : sampling frequency [Hz]
 --- Optional parameter for EMG onset
 parm.p_val : P-value corresponding to the threshold   [0.001]
 parm.t_smooth : moving average window length          [20 ms]

 --- Output
 y0 = z0 : threshold value estimated from moving average
 z  : smoothed EMG power by TKE operator

 2009-09-05 Masa-aki Sato
 2012-1-8   Masa-aki Sato  set z0 = y0
 2012-2-18  Masa-aki Sato  use TKE operator

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function    [y0, z0, z] = vb_get_emg_threshold(y,parm)
0002 % Estimate threshold for smoothed EMG power
0003 %    [y0, z0, z] = vb_get_emg_threshold(y,parm)
0004 %
0005 % --- Input
0006 % y             : EMG signal
0007 %
0008 % parm.fsamp    : sampling frequency [Hz]
0009 % --- Optional parameter for EMG onset
0010 % parm.p_val : P-value corresponding to the threshold   [0.001]
0011 % parm.t_smooth : moving average window length          [20 ms]
0012 %
0013 % --- Output
0014 % y0 = z0 : threshold value estimated from moving average
0015 % z  : smoothed EMG power by TKE operator
0016 %
0017 % 2009-09-05 Masa-aki Sato
0018 % 2012-1-8   Masa-aki Sato  set z0 = y0
0019 % 2012-2-18  Masa-aki Sato  use TKE operator
0020 %
0021 % Copyright (C) 2011, ATR All Rights Reserved.
0022 % License : New BSD License(see VBMEG_LICENSE.txt)
0023 
0024 % Default value
0025 mode = 1;  % low-pass filtering
0026 t_smooth  = 20;  % 10/20 ms average (100/50 Hz)
0027 p_val = [0.001];
0028 
0029 if isfield(parm,'t_smooth'), t_smooth = parm.t_smooth;end
0030 if isfield(parm,'p_val'),    p_val  = parm.p_val ;end
0031 
0032 % sampling frequency
0033 fs  = parm.fsamp ; 
0034 
0035 switch    mode
0036 case    0
0037     % smoothed signal for threshold estimation
0038     tau  = round(t_smooth * fs/ 1000);  % msec -> sample number
0039     yave = filter( ones(tau,1)/tau, 1, abs(y));
0040 case    1
0041     % low-pass filtering
0042     Norder = 2;
0043     [B,A] = butter(Norder, 1000/(t_smooth * fs /2) ,'low');
0044     yave  = filtfilt(B,A, abs(y) );
0045 end
0046 
0047 % Estimate threshold by gamma distribution approximation
0048 [y0, A ,hy ,ylist] = vb_gamma_dist_param(yave, p_val(1));
0049 
0050 % TKE operator EMG filter
0051 z = vb_emg_filter(y, parm.fsamp);
0052 
0053 % set (smoothed EMG threshold) = (EMG threshold)
0054 z0 = y0;
0055 
0056 if isfield(parm,'plot_mode') && parm.plot_mode > 0
0057     vb_show_hist_threshold(y0, A ,hy ,ylist);
0058     Nbin = 1000;
0059     [hz, zx] = hist(z,Nbin);
0060     figure;
0061     plot(zx,hz);
0062     hold on
0063     plot([z0 z0],[0 max(hz)*0.2],'r-')
0064     xlim([0 z0*3]);
0065 end
0066 
0067 return
0068 % ---- END of program ----