% function activationTimeMap(A)
close all
[rows,cols,frames] = size(A);
nPix = rows*cols;
%B = zeros(nPix, frames);
% for i = 1:frames
%     a = A(:,:,i);
%     B(:,i) = a(:);
% end
B = reshape(A,[nPix,frames]);
Fs = 50;
%%
X = fft(B');



P2 = abs(X/frames);
P1 = P2(1:frames/2+1,:);
P1(2:end-1,:) = 2*P1(2:end-1,:);

f = Fs*(0:(frames/2))/frames;
% plot(f,P1) 
% title('Single-Sided Amplitude Spectrum of X(t)')
% xlabel('f (Hz)')
% ylabel('|P1(f)|')
%%
Pavg = mean(P1,2);
Pavg(1) = 0; % no DC component

[~,ind] = max(Pavg);
% fmax = f(ind);
%% 
%for ind = 1:108
fftMaxFreq = X(ind,:);
fftAbs = abs(fftMaxFreq);
fftAngle = angle(fftMaxFreq);

%colormap(flipud(parula))
C = reshape(fftAngle(:),[rows,cols]);

E = reshape(fftAbs(:),[rows,cols]);
% figure
% imagesc(real(C)) 
% figure
% imagesc(imag(C))
% midp = 64;
% C = medfilt2(C);
% CL = unwrap(C(:,midp:-1:1),[],2);
% CR = unwrap(C(:,midp+1:end),[],2);
% C = [fliplr(CL),CR];
figure(1)

C = medfilt2(C, [5 5],'symmetric');
C = C/(2*pi*f(ind));
C = C - max(C(:));

C_C = imgaussfilt(C,50);
C_C = C_C - min(C_C(:));
hIm = imagesc(C,[-.5, 0]); hold on
colorbar
contour(C_C, 0:.05:1, 'r','LineWidth',2); hold off
% s1 = sprintf('phase of signal at f=%4.2fHz (radians)',f(ind));
s1 = sprintf('time delay at f=%4.2fHz (s)',f(ind));
title(s1)

% figure(2)
% imagesc(E)
% s2 = sprintf('intensity of signal at f=%4.2fHz',f(ind));
% title(s2)
% colorbar
% pause 
% end  
