% decomp - Temporally decompose single-channel data, such as 
%          electroencephalography (EEG) data, during reaction time 
%          tasks into a stimulus- and response-locked components using 
%          reaction times.
%
% Usage: 
%         >> [sc,rc] = decomp(data,rt,stim);
%
% Input:
%  data - Single-channel data to decompose.
%         Format (frames,trials)
%  rt   - Number of data points between stimulus and response onset.
%         Format (trials,1)
%  stim - Number of data points before stimulus onset.
%
% Output:
%  sc  - Extracted stimulus-locked component.
%         Format (frames,1)
%         Number of data points before stimulus onset is 'stim' in 
%         the input.
%  rc  - Extracted response-locked component.
%         Format (frames,1)
%         Number of data points before response onset is 'stim' in 
%         the input.
%
% Author: Yusuke Takeda
% Educational Physiology Laboratory
% Graduate School of Education, The University of Tokyo
% 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan
%
% Copyright (C) 2007 Yusuke Takeda (takeda@p.u-tokyo.ac.jp)
%
% Last modified: Oct 11, 2007
%
% This program is free software; you can redistribute it and/or 
% modify it under the terms of the GNU General Public License as 
% published by the Free Software Foundation; either version 2 of the 
% License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%

function [sc,rc]=decomp(data,rt,stim);

%%%%%%%%%%%%%% Check inputs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin < 3;
  error('You should provide three inputs')
end

[Frames,Trials] = size(data); 
if min(rt) < 1 | max(rt) > Frames;
  error('rt should be between 1 and the size of the first dimension of data');
end

if length(rt) ~= Trials;
 error('The length of rt should be the same as the size of the second dimension of data');
end;

if stim < 1 | stim > Frames;
 error('stim should be between 1 and the size of the first dimension of data');
end

%%%%%%%%%%%%%% Initialise parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c = [1;zeros(Frames-1,1)]; 
rt1 = zeros(Trials,1);
rt2 = zeros(Frames,1);
data1 = zeros(Frames,Trials);
data2 = zeros(Frames,Trials);

%%%%%%%%%%%%%% Sort trials by RTs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rts=sortrows([(1:Trials)',rt],2);
for I = 1:Trials;
  rt1(I,1) = rts(I,2);
  data1(:,I) = data(:,rts(I,1));
end;

%%%%%%%%%%%%%% Take a discrete Fourier transform %%%%%%%%%%%%%%%%%%%%
j = 1;
k = 1;
while k <= Trials;
  a = find(rt1==rt1(k,1));
  n_a = size(a,1);
  data2(:,j) = sum(data1(:,a(1,1):a(n_a,1)),2);
  rt2(rt1(k,1),j) = n_a;
  k = a(n_a,1)+1;
  j = j+1;
end;
ns = j-1;
data2 = fft(data2(:,1:ns));
rt2 = fft(rt2(:,1:ns));

%%%%%%%%%%%%%% Average data across trials %%%%%%%%%%%%%%%%%%%%%%%%%%%
data2_ave = sum(data2,2)/Trials;
rt2_ave = sum(rt2,2)/Trials;

%%%%%%%%%%%%%% Solve equations simultaneously %%%%%%%%%%%%%%%%%%%%%%%
er = zeros(Frames,1);
es = zeros(Frames,1);
for I = 1:ns;
  cc = rt2(:,I)/rt2(1,I)-rt2_ave+c;
  es = (rt2(:,I).*data2_ave-rt2_ave.*data2(:,I))./cc+es;
  er = (data2(:,I)-data2_ave*rt2(1,I))./cc+er;
end;

%%%%%%%%%%%%%% Take a inverse discrete Fourier transform %%%%%%%%%%%%
sc = real(ifft(es/Trials));
rc = real(ifft(er/Trials));

%%%%%%%%%%%%%% End %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
return
