function [Y] = meeg_decode_circ(X,method)

% [Y] = meeg_decode_circ(X,method)
%
% X      - cell vector (of conditions) each containing a trial x time matrix of phase angles
% method - 'SI' (use synchronization index) or 'DIST' (use circular distance, default)
%
% Y - decoding/confusion matrix (proportion), with columns refering to the 
%     actual category/conditions and rows to the decoded ones
%
% Reference:
% Cohen MX (2008) Assessing transient cross-frequency coupling in EEG data.
%		Journal of Neuroscience Methods 168:494-499.
% Ng BS, Logothetis NK, Kayser C (in press) EEG Phase Patterns Reflect the Selectivity of Neural Firing.
%		Cerebral Cortex.
%
% Description: The script calculates a decoding/confusion matrix for
% several circular conditions. In particular, a template for each condition
% is calculated (circular mean), and then each single trial is compared to
% the templates and classified as belonging to a particular condition (SI or DIST).
% The single trial is excluded from the template in the corresponding instance.
%
% ---------
%
%    Copyright (C) 2012, B. Herrmann
%    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 3 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.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
% ------------------------------------------------------------------------------
% B. Herrmann, Email: herrmann.b@gmail.com, 2012-11-06

if nargin < 1 || isempty(X), error('X needs to be provided!'); end
if nargin < 2 || isempty(method), method = 'DIST'; end


Y = zeros(length(X),length(X));
for ii = 1 : length(X)
	% get circular mean across trials for all conditions (i.e. establish templates)
	T = zeros(length(X),size(X{1},2));
	for kk = 1 : length(X)
		T(kk,:) = angle(sum(exp(1i*X{kk}),1));
	end
	
	% do classification for each trial
	for jj = 1 : size(X{ii},1)
		% get single trial and adjusted template
		x = X{ii}(1:size(X{ii},1)==jj,:);
		T(ii,:) = angle(sum(exp(1i*X{ii}(1:size(X{ii},1)~=jj,:))));
		
		if strcmp(method,'SI')
			% get synchronization index
			SI = abs(mean(exp(1i*bsxfun(@minus,T.',x.'))));

			% get index of largest synchronization
			[~, ix] = max(SI);
		elseif strcmp(method,'DIST')
			% get circular distance between templates and trial
			dist = angle(exp(1i*T) ./ repmat(exp(1i*x),size(T,1),1));

			% get circular mean across time and then the smallest difference from zero
			[~, ix] = min(abs(angle(sum(exp(1i*dist),2))));
		else
			error('Unrecognized method!!!')
		end
		% add to decoding / confusion matrix
		Y(:,ii) = Y(:,ii) + (1:size(T,1)==ix)';
	end
	% get decoding / confusion matrix in proportion
	Y(:,ii) = Y(:,ii) / size(X{ii},1);
end


