clear all;


%% Create a complex wavelet

% Params:
freq=5;                  % frequency of the sine wave (Hz)
sampRate=512;            % sampling rate

% Centered time vector
time = -2:1/sampRate:2;
% make a wavelet
c = 5;
s = c/(2*pi*freq);
cmw = exp( 1i*2*pi*freq(*time) .* exp(-time.^2./(2*s^2));

%% TASK

% Params:
freq=5;                  % frequency of the sine wave (Hz)
sampRate=512;            % sampling rate

% Centered time vector
time = -2:1/sampRate:2;
% make a wavelet
c = linspace(1,10,11);

for ci = 1:size(c,2)
    s = c(ci)/(2*pi*freq);
    cmw(ci,:) = exp( 1i*2*pi*freq*time) .* exp(-time.^2./(2*s^2));
end

freq = 2:1:10
for ci = 1:length(freq)
    s = 5/(2*pi*freq(ci));
    cmw(ci,:) = exp( 1i*2*pi*freq(ci)*time) .* exp(-time.^2./(2*s^2));
end

figure(4),clf
for ci = 1:length(c)
    subplot(3,4,ci)
        plot(time,real(cmw(ci,:)),LineWidth=3)
end

% 3-D
plot3(time,real(cmw(3,:)),imag(cmw(3,:)),LineWidth=3)
%%

load("FCz_sample.mat")

% Convolution and FFT parameters:
kernel = length(time);
htime = (length(time)-1)/2;
ndata  = length(signal);
nconv  = kernel+ndata-1;

% FFTs of kernel and signal
cx = fft(cmw,nconv);
dx = fft(signal,nconv);
% Normalize kernel (optional)
cx = cx./max(cx);

% Convolution
rc = cx .* dx;
as = ifft(rc,nconv);
% % Cut off wings
as = as(htime+1:end-htime);
% % Take a magnitude of the complex signal
as = abs(as).^2;

%% T2

% Load the sample data:
load("FCz_sample.mat")

% Wavelet parameters
sampRate = 512;       
mtime = -2:1/sampRate:2;
htime = (length(mtime)-1)/2;
num_frex    = 30;
min_freq    = 1;
max_freq    = 30;
frex        = logspace(log10(min_freq),log10(max_freq),num_frex); 
rcycles     = [3 10];
ncycles     = logspace(log10(rcycles(1)),log10(rcycles(end)),num_frex);


% Convolution and FFT parameters
kernel  = length(mtime);
ndata   = length(signal);
nconv   = kernel+ndata-1;
timevec = linspace(0,(ndata-3)/sampRate,ndata) * 1000;
% Create a set of wavelets
cmwFFT  = zeros(length(frex),nconv);
for fi=1:num_frex
    s = ncycles(fi)/(2*pi*frex(fi));
    cmw = exp(2*1i*pi*frex(fi).*mtime) .* exp(-mtime.^2./(2*s^2)); %complex morlet wavelet creation

    tmp = fft(cmw,nconv);
    cmwFFT(fi,:) = tmp./max(tmp); %uV normalization units  
end

dataFFT     = fft(signal,nconv);
baseidx     = dsearchn(timevec',[1000 2000]');
[tf_result,bs_result] = deal( zeros(size(frex,2),ndata) );

for j = 1:length(frex)
    % Convolution across frequencies
    as = ifft(cmwFFT(j,:).*dataFFT,nconv);
    as = as(htime+1:end-htime); 
    % Extract power values
    tf_result(j,:) = abs(as).^2; 
    % Baseline normalization
    tf_baseline = mean(tf_result(j,baseidx(1):baseidx(2)),2);    
    bs_result(j,:) = 100*( (tf_result(j,:)- tf_baseline )./tf_baseline);
end

% Plot the result
figure; clf;set(gcf, 'color','w', 'WindowState', 'maximize');
contourf(timevec,frex,bs_result,200,'linecolor','none','fill','on'),hold on 
set(gca,'ylim',ylim,'xlim',xlim,'clim',[-1000 1000],'FontSize',35,...
    'LineWidth',3) 
    box off;colormap (jet);colorbar