Home > vbmeg > functions > device > trigger_timing > vb_emg_spectral_check.m

vb_emg_spectral_check

PURPOSE ^

Check EMG frequency spectrum

SYNOPSIS ^

function [trig, i_emg, sp, f] = vb_emg_spectral_check(y,trig,parm)

DESCRIPTION ^

 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)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Mon 22-May-2023 06:53:56 by m2html © 2005