clear all;
format compact

modrate = 2:0.5:7;
dur     = 3;
Sf      = 200;
nTrials = 100;

rand('seed',645)

[signals signalsT] = deal({});
t = 0:1/Sf:dur-1/Sf;
for ii = 1 : nTrials
	for jj = 1 : length(modrate)
		phi = (rand(1)-0.5)*pi/1.5;
		tmpsignal = sin(pi*2*modrate(jj)*t+phi);
		noise     = ((rand([1 length(tmpsignal)])-0.5)*2)*5;
		signals{jj}(ii,:) = angle(hilbert(tmpsignal + noise));
		signalsT{jj}(ii,:) = tmpsignal + noise;
	end
end

idxplot = [1 4];
% plot time-domain
figure, set(gcf,'Color',[1 1 1])
nn = 1;
for ii = idxplot
	subplot(2,2,nn)
	plot(t, mean(signalsT{ii}),'-k')
	xlim([0 dur])
	nn = nn + 1;
	title(['Mod Rate: ' num2str(modrate(ii)) 'Hz'])
end
for ii = idxplot
	subplot(2,2,nn)
	plot(t, signalsT{ii}(1,:),'-k')
	xlim([0 dur])
	nn = nn + 1;
end

% plot phase time series
figure, set(gcf,'Color',[1 1 1])
nn = 1;
for ii = idxplot
	subplot(2,2,nn)
	plot(t, circ_mean(signals{ii}),'-k')
	xlim([0 dur])
	nn = nn + 1;
	title(['Mod Rate: ' num2str(modrate(ii)) 'Hz'])
end
for ii = idxplot
	subplot(2,2,nn)
	plot(t, signals{ii}(1,:),'-k')
	xlim([0 dur])
	nn = nn + 1;
end



% do classification
Y = meeg_decode_circ(signals,'DIST');

% plot classification results
figure, set(gcf,'Color',[1 1 1])
set(gca,'FontSize',12)
imagesc(Y)
set(gca,'XTick',1:length(modrate),'XTickLabel',modrate)
set(gca,'YTick',1:length(modrate),'YTickLabel',modrate)
set(gca,'CLim',[0 0.4])
daspect([1 1 1])
colormap(flipud(pink(201)))
colorbar
xlabel('Decoded modulation rate')
ylabel('Stimulus modulation rate')

CC = diag(Y);
figure, set(gcf,'Color',[1 1 1]), hold on
set(gca,'FontSize',12)
for ii = 1 : length(CC)
	rectangle('Position',[ii-0.45 0 0.9 CC(ii)],'FaceColor',[0.7454 0.4774 0.4714])
end
plot([0 12],[1/length(CC) 1/length(CC)],':k')
ylim([0 0.4])
set(gca,'XTick',1:length(modrate),'XTickLabel',modrate)
xlabel('Modulation rate')
ylabel('Decoding accuracy')

