Home > bioelectromagnetism > eeg_amplitudearea.m

eeg_amplitudearea

PURPOSE ^

eeg_amplitudearea() - Resamples an ERP average using spline interpolation

SYNOPSIS ^

function [channels,overall_amplitude] = eeg_amplitudearea(EEG, channels, resrate, wstart, wend)

DESCRIPTION ^

 eeg_amplitudearea() - Resamples an ERP average using spline interpolation 
                       at a new sample rate (resrate) in Hz to get the exact limits 
                       of the window of integration. Finely samples the window 
                       and adds together very narrow rectangles capped by 
                       right-angled triangles under the window. Output is in uV. 
                       Trade-off between speed and number of resamples and number of 
                        channels selected occurs.
 Usage:
      >> [channels, amplitude] = eeg_amplitudearea(EEG,channels, resrate, wstart, wend);
 Inputs:
   EEG         - EEGLAB data struct containing a (3-D) epoched data matrix 
   channels    - vector of channel indices
   resrate     - resampling rate for window of integration in Hz
   wstart      - start of window of integration in ms post stimulus-onset
   wend        - end of window of integration in ms post stimulus-onset

 Outputs:
   channels    - a vector of channel indices.
   amplitude   - 1-dimensional array in uV for the channels

 Example
    >> [channels, amplitude] = eeg_amplitudearea(EEG,[12 18 25 29], 2000, 90.52, 120.52);

 Author: Tom Campbell, Helsinki Collegium for Advanced Studies, Biomag Laboratory, 
         Engineering Centre, Helsinki University Central Hospital Helsinki Brain 
         Research Centre (tom.campbell@helsinki.fi) Spartam nanctus es: Hanc exorna. 
         Combined with amplitudearea_msuV() by Darren Weber, UCSF 28/1/05
         Retested and debugged Tom Campbell 2/2/05

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % eeg_amplitudearea() - Resamples an ERP average using spline interpolation
0002 %                       at a new sample rate (resrate) in Hz to get the exact limits
0003 %                       of the window of integration. Finely samples the window
0004 %                       and adds together very narrow rectangles capped by
0005 %                       right-angled triangles under the window. Output is in uV.
0006 %                       Trade-off between speed and number of resamples and number of
0007 %                        channels selected occurs.
0008 % Usage:
0009 %      >> [channels, amplitude] = eeg_amplitudearea(EEG,channels, resrate, wstart, wend);
0010 % Inputs:
0011 %   EEG         - EEGLAB data struct containing a (3-D) epoched data matrix
0012 %   channels    - vector of channel indices
0013 %   resrate     - resampling rate for window of integration in Hz
0014 %   wstart      - start of window of integration in ms post stimulus-onset
0015 %   wend        - end of window of integration in ms post stimulus-onset
0016 %
0017 % Outputs:
0018 %   channels    - a vector of channel indices.
0019 %   amplitude   - 1-dimensional array in uV for the channels
0020 %
0021 % Example
0022 %    >> [channels, amplitude] = eeg_amplitudearea(EEG,[12 18 25 29], 2000, 90.52, 120.52);
0023 %
0024 % Author: Tom Campbell, Helsinki Collegium for Advanced Studies, Biomag Laboratory,
0025 %         Engineering Centre, Helsinki University Central Hospital Helsinki Brain
0026 %         Research Centre (tom.campbell@helsinki.fi) Spartam nanctus es: Hanc exorna.
0027 %         Combined with amplitudearea_msuV() by Darren Weber, UCSF 28/1/05
0028 %         Retested and debugged Tom Campbell 2/2/05
0029 
0030 
0031 function [channels,overall_amplitude] = eeg_amplitudearea(EEG, channels, resrate, wstart, wend)
0032 
0033 if wstart > wend
0034     error ('ERROR: wstart must be greater than wend')
0035 else
0036     [channels, overall_amplitude] = amplitudearea_msuV (EEG,channels, resrate, wstart, wend)
0037     overall_amplitude = overall_amplitude/(wend - wstart)
0038 end
0039 
0040 return
0041 
0042 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0043 function [channels, overall_area] = amplitudearea_msuV (EEG, channels, resrate, wstart, wend)
0044 
0045 %if ndim(EEG.data) ~= 3
0046 %  error('EEG.data must be 3-D data epochs');
0047 %end
0048 erp = mean(EEG.data,3);
0049 
0050 [tmp ind1] =min( abs( EEG.times - wstart ) ); % closest time index to wstart
0051 [tmp ind2] =min( abs( EEG.times - wend ) ); % closest time index to wend
0052 restep = 1000/resrate
0053 
0054 
0055 if EEG.times(ind1) > wstart 
0056     ind1= ind1 -1;
0057 end
0058 
0059 if EEG.times(ind2) < wend 
0060     ind2= ind2 +1;
0061 end    
0062 
0063 for x= ind1:ind2
0064     t = (x -ind1)+1;
0065     tim(t) = EEG.times(x);
0066 end
0067 
0068 tr = 1
0069 timr(tr) = wstart;
0070 while timr(tr) < wend
0071     tr = tr + 1;
0072     timr(tr) = timr(tr-1)+ restep;
0073 end
0074 
0075 for x = 1:size(channels,2)
0076     channel = channels(x)
0077     %resamples
0078     rerp(x, 1:tr) = spline(tim(:),erp(channel, ind1:ind2), timr(1:tr))
0079     for y = 1:(tr -1)
0080         %identify which sample is the height of the rectangle under the
0081         %curve and which sample is the height of the triangle capping the
0082         %rectangle
0083         if (abs(rerp(x,y)) <= abs(rerp(x,y+1)))
0084             rhtsamp = y
0085             thtsamp = y + 1
0086         else 
0087             rhtsamp = y + 1
0088             thtsamp = y
0089         end
0090         if y < (tr-1) | ((y == (tr-1)) & (wend - timr(tr - 1)) == restep) 
0091             if ((rerp(x,y) > 0)& (rerp(x,y+1) < 0))|((rerp(x,y) < 0)& (rerp(x,y+1) > 0))
0092                 %handles samples containing a zerocrossing with two
0093                 %triangles and trigonometry
0094                 opposite = abs(rerp(x,y))+abs(rerp(x,y+1))
0095                 adjacent = restep
0096                 tantheta = opposite/adjacent
0097                 zerocross(x,y) = abs(rerp(x,y))/tantheta
0098                 remainder(x,y) = adjacent - zerocross(x,y)
0099                 area(x,y) = ((zerocross(x,y)/2)* (rerp(x,y)))+ ((remainder(x,y)/2)* (rerp(x,y+1)))
0100             else       
0101                 rectanglearea = rerp(x,(rhtsamp))*restep 
0102                 trianglearea = (rerp(x,(thtsamp))-rerp(x,(rhtsamp)))*restep/2 
0103                 area(x,y) =  rectanglearea  + trianglearea
0104             end 
0105         elseif (y == (tr-1)) & ((wend - timr(tr - 1)) < restep)
0106             endstep = wend - timr(tr - 1)
0107             opposite = abs(rerp(x,y))+abs(rerp(x,y+1))
0108             adjacent = restep
0109             tantheta = opposite/adjacent
0110             zerocross = abs(rerp(x,y))/tantheta
0111              if (((rerp(x,y) > 0)& (rerp(x,y+1) < 0))|((rerp(x,y) < 0)& (rerp(x,y+1) > 0)))&(zerocross < endstep)
0112                 %handles samples containing a zerocrossing with two
0113                 %triangles and trigonometry
0114                 remainder = endstep - zerocross
0115                 opposite2 = tantheta * remainder * sign(rerp(x,y)) * -1
0116                 area(x,y) = ((zerocross/2)* (rerp(x,y)))+ ((remainder/2)* opposite2)
0117             elseif thtsamp > rhtsamp
0118                 % test for amplitude increase from left to right and fits
0119                 % area with a rectangle capped by a triangle
0120                 opposite = abs(rerp(x,(thtsamp))) - abs(rerp(x,(rhtsamp)))
0121                 adjacent = restep
0122                 tantheta = opposite/adjacent
0123                 triht = tantheta * endstep * sign(rerp(x,(thtsamp)))
0124                 trianglearea = (endstep/2)*triht 
0125                 rectanglearea = rerp(x,(rhtsamp))*endstep
0126                 area(x,y) =  rectanglearea + trianglearea
0127             elseif rhtsamp > thtsamp 
0128                  % tests for amplitude decrease from left to right and fits area with
0129                 % two rectangles capped by a triangle with trigonometry
0130                 opposite = abs(rerp(x,(thtsamp))) - abs(rerp(x,(rhtsamp))) 
0131                 adjacent = restep
0132                 tantheta = opposite/adjacent
0133                 excessstep = restep - endstep
0134                 extrarectangleht = tantheta* excessstep * sign(rerp(x,(thtsamp)))
0135                 extrarectanglearea = extrarectangleht * endstep *sign(rerp(x,(thtsamp)))
0136                 triht = (opposite - extrarectangleht) * sign(rerp(x,(thtsamp)))
0137                 trianglearea = (endstep/2)*triht *sign(rerp(x,(thtsamp)))
0138                 rectanglearea = rerp(x,(rhtsamp))*endstep
0139                 area(x,y) =  rectanglearea + extrarectanglearea + trianglearea
0140             end
0141         end
0142     end
0143 end
0144 
0145 for x = 1:(size(channels,2))
0146     overall_area(x) = sum(area(x,:))
0147 end
0148 
0149 return
Generated on Mon 15-Aug-2005 15:36:19 by m2html © 2003