vb_emg_spectral_check

 Check EMG frequency spectrum
 If power of EMG frequency range (above 10 Hz) is 
 less than lower frequency power, it is rejected

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

0001 function    [trig, i_emg, sp, f] = vb_emg_spectral_check(y,trig,parm)
0002 % Check EMG frequency spectrum
0003 % If power of EMG frequency range (above 10 Hz) is
0004 % less than lower frequency power, it is rejected
0005 %
0006 % Copyright (C) 2011, ATR All Rights Reserved.
0007 % License : New BSD License(see VBMEG_LICENSE.txt)
0008 
0009 T1 = 50;  % Pretrigger   100 ms
0010 T2 = 200; % Posttrigger  200 ms
0011 
0012 Fhigh = 100; % upper frequency [Hz]
0013 Flow  = 10;  % lower threshold frequency [Hz]
0014 
0015 % Threshold value for low/high ratio
0016 if isfield(parm,'emg_low_ratio') && ~isempty(parm.emg_low_ratio)
0017     R0 = parm.emg_low_ratio;
0018 else
0019     R0 = 1;
0020 end
0021 
0022 N  = length(trig);
0023 T0 = length(y);
0024 
0025 fs = parm.fsamp;
0026 tw = ceil( [T1 T2]* fs/1000);
0027 
0028 % Time index for one trial
0029 t  = [-tw(1):tw(2)-1];
0030 T  = length(t);
0031 
0032 % Time index for all trials extracted from 'trig'
0033 t_trial = repmat(t(:),[1 N]) + repmat(round(trig(:))',[T 1]);
0034 t_trial = max(t_trial,1);
0035 t_trial = min(t_trial,T0);
0036 
0037 % signal for all trials
0038 xdata = y(t_trial);
0039 
0040 % Frequency list for FFT
0041 fd = fs/T;
0042 f  = 0:fd:(fs/2);
0043 
0044 % Window function
0045 xx  = vb_repmultiply(xdata ,hanning(T));
0046 
0047 % Frequency spectrum
0048 sp = abs(fft(xx)).^2;
0049 
0050 % remove hifh frequency than 'Fhigh'
0051 jx = find( f <= Fhigh);
0052 f  = f(jx);
0053 sp = sp(jx,:);
0054 
0055 % Lower frequency component than 'Flow'
0056 lx = find( f <= Flow );
0057 
0058 % ratio of low/high component
0059 rs = sum(sp(lx,:),1)./sum(sp((lx(end)+1):end,:),1);
0060 
0061 i_emg = find( rs <  R0); % index for correct EMG
0062 
0063 trig = trig(i_emg);
0064 
0065 i_rej = find( rs >= R0); % index for rejection
0066 
0067 if ~isempty(i_rej)
0068     figure;
0069     subplot 211
0070     plot(xdata(:,i_rej))
0071     title('Rejected EMG')
0072     xlabel('Time')
0073     ylabel('EMG')
0074     
0075     subplot 212
0076     plot(f,sp(:,i_rej))
0077     title('Rejected EMG spectrum')
0078     xlabel('Frequency [Hz]')
0079     ylabel('Spectral power')
0080 end