eeg_amplitudearea2

 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
         Reconceived, factored somewhat, tested and debugged Tom Campbell 13:24 23.3.2005

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 %         Reconceived, factored somewhat, tested and debugged Tom Campbell 13:24 23.3.2005
0030 
0031 
0032 
0033 function [channels,overall_amplitude] = eeg_amplitudearea2(EEG, channels, resrate, wstart, wend)
0034 
0035 if wstart > wend
0036     error ('ERROR: wstart must be greater than wend')
0037 else
0038     [channels, overall_amplitude] = eeg_amplitudearea_msuV (EEG,channels, resrate, wstart, wend)
0039     overall_amplitude = overall_amplitude/(wend - wstart)
0040 end
0041 
0042 return
0043 
0044 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0045 function [channels, overall_area] = eeg_amplitudearea_msuV (EEG, channels, resrate, wstart, wend)
0046 
0047 %if ndim(EEG.data) ~= 3
0048 %  error('EEG.data must be 3-D data epochs');
0049 %end
0050 erp = mean(EEG.data,3);
0051 
0052 [tmp ind1] =min( abs( EEG.times - wstart ) ); % closest time index to wstart
0053 [tmp ind2] =min( abs( EEG.times - wend ) ); % closest time index to wend
0054 restep = 1/resrate
0055 
0056 
0057 if EEG.times(ind1) > wstart 
0058     ind1= ind1 -1;
0059 end
0060 
0061 if EEG.times(ind2) < wend 
0062     ind2= ind2 +1;
0063 end    
0064 
0065 for x= ind1:ind2
0066     t = (x -ind1)+1;
0067     tim(t) = EEG.times(x);
0068 end
0069 
0070 tr = 1
0071 timr(tr) = wstart;
0072 while timr(tr) < wend
0073     tr = tr + 1;
0074     timr(tr) = timr(tr-1)+ restep;
0075 end
0076 
0077 
0078 for x = 1:size(channels,2)
0079     channel = channels(x)
0080     %resamples
0081     rerp(x, 1:tr) = spline(tim(:),erp(channel, ind1:ind2), timr(1:tr))
0082     pent = timr(tr - 1)
0083     overall_area(x) = 0
0084     for y = 1:(tr -1)
0085         v1 =  rerp(x,(y))
0086         v2 =  rerp(x,(y+1))
0087         if ((v1 > 0) & (v2 < 0)) | ((v1 < 0) & (v2 > 0))
0088             if (y == (tr-1)) & (timr(y+1)> wend)
0089                 area1 = zero_crossing_truncated(v1, v2, restep, wend, pent)    
0090             else    
0091                 area1 = zero_crossing(v1, v2, restep)
0092             end
0093         else
0094             if( y == (tr-1)) & (timr(y+1)> wend)
0095                 area1 = rect_tri_truncated(v1, v2, restep,wend,pent)
0096             else
0097                 area1 = rect_tri(v1, v2, restep)
0098             end
0099         end
0100         overall_area(x) = overall_area(x) + area1
0101     end
0102 end
0103 
0104 return
0105 
0106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0107 function [area] = zero_crossing(v1,v2,step)
0108     if (v1 > v2)
0109         T1 = v1
0110         T2 = v2
0111     else
0112         T1 = v2
0113         T2 = v1
0114     end
0115     tantheta = (abs(T1)+ abs(T2))/step
0116     if (v1 > v2)
0117         %decline
0118         z = abs(T1)/tantheta
0119         tr1= abs(T1)*(z/2)
0120         tr2= abs(T2)*((step-z)/2)
0121     else
0122         %incline
0123         z = abs(T2)/tantheta
0124         tr2= abs(T2)*(z/2)
0125         tr1= abs(T1)*((step-z)/2)
0126     end
0127     [area] = (tr1 - tr2)
0128 return
0129 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0130 function [area] = zero_crossing_truncated(v1,v2,step,wend,pent)
0131     if (v1 > v2)
0132         T1 = v1
0133         T2 = v2
0134     else
0135         T1 = v2
0136         T2 = v1
0137     end
0138     tantheta = (abs(T1)+ abs(T2))/step
0139     s = wend - pent
0140     if (v1 > v2)
0141         z = abs(T1)/tantheta
0142         if s < z   
0143             %decline,truncated before zerocrossing
0144             t1 = tantheta * s
0145             r1 = abs(T1)-abs(t1)
0146             tr1= abs(t1)*(s/2)
0147             tr2= 0
0148             rect1 = r1*s
0149             rect2 = 0
0150         else
0151             %decline,truncated after zerocrossing
0152             t2= tantheta*(s-z)
0153             tr1= abs(T1)*(z/2)
0154             tr2 = abs(t2)*((s-z)/2)
0155             rect1 = 0
0156             rect2 = 0
0157         end    
0158     else
0159         z = abs(T2)/tantheta
0160         if s < z
0161             %incline,truncated before zerocrossing
0162             t2 = tantheta * s
0163             r2 = abs(T2)-abs(t2)
0164             tr1= 0
0165             tr2= abs(t2)*(s/2)            
0166             rect1 = 0
0167             rect2 = r2*s
0168         else
0169             %incline,truncated after zerocrossing
0170             t1= tantheta*(s-z)
0171             tr1 = abs(t1)*((s-z)/2)
0172             tr2 = abs(T2) * (z/2)
0173             rect1 = 0
0174             rect2 = 0
0175         end
0176     end
0177 
0178 [area] = ((rect1 + tr1) - (rect2 + tr2))
0179 return
0180 
0181 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0182 function [area] = rect_tri(v1,v2,step)
0183     if (abs(v1) > abs(v2))
0184         T = abs(v1)-abs(v2)
0185         R = abs(v2)
0186     else
0187         T = abs(v2)-abs(v1)
0188         R = abs(v1)
0189     end
0190     rect = R*step
0191     tri = T*(step/2)
0192     if v1 > 0
0193         area = 1* (rect+tri)
0194     else
0195         area = -1 * (rect+tri)
0196     end
0197 return
0198 
0199 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0200 function [area] = rect_tri_truncated(v1,v2,step,wend,pent)
0201     if (abs(v1) > abs(v2))
0202         T = abs(v1)-abs(v2)
0203         R = abs(v2)
0204     else
0205         T = abs(v2)-abs(v1)
0206         R = abs(v1)
0207     end
0208     
0209     tantheta = abs(T)/step
0210     s = wend -pent
0211     
0212     if (v1>0)
0213         if v1 >v2
0214         %positive decline
0215             t = tantheta*s
0216             e = abs(T)-abs(t)
0217             rect = s*R
0218             exrect = s*e
0219             tri = (s/2)*R
0220         else
0221         %positive incline
0222             t = tantheta*s
0223             rect = s*R
0224             exrect = 0
0225             tri = (s/2)*R
0226         end
0227     else
0228        if v1 >v2
0229         %negative decline
0230             t = tantheta*s
0231             rect = s*R
0232             exrect = 0
0233             tri = (s/2)*R
0234         else
0235         %negative incline
0236             t = tantheta*s
0237             e = abs(T)-abs(t)
0238             rect = s*R
0239             exrect = s*e
0240             tri = (s/2)*R
0241         end
0242     end
0243     tri = T*(step/2)
0244     if v1 > 0
0245         area = 1* (rect+exrect+tri)
0246     else
0247         area = -1 * (rect+exrect+tri)
0248     end
0249 return