%% Power analysis script

% Shows the power spectrum for each mouse
% - script is used mainly to confirm that the algorithms are working as
% intended, mainly showing that there ocurrs circadian rhythm
% - no further spectral analyses were carried out


close all;
clear;

%% LOAD

load data\travelled\dlc_travelled_jc20181211_1_05.mat
dlc_travelled_jc20181211_1_05 = DLC_travelled;
travelled_ourDLC_5fps = medfilt1(dlc_travelled_jc20181211_1_05.travelled_ourDLC, 10*5)*5/100;
dlc_travelled_jc20181211_1_05.travelled_ourDLC = travelled_ourDLC_5fps;

load data\travelled\dlc_travelled_jc20181218_1_05.mat
dlc_travelled_jc20181218_1_05 = DLC_travelled;
travelled_ourDLC_5fps = medfilt1(dlc_travelled_jc20181218_1_05.travelled_ourDLC, 10*5)*5/100;
dlc_travelled_jc20181218_1_05.travelled_ourDLC = travelled_ourDLC_5fps;

load data\travelled\dlc_travelled_jc20181219_2_05.mat
dlc_travelled_jc20181219_2_05 = DLC_travelled;
travelled_ourDLC_5fps = medfilt1(dlc_travelled_jc20181219_2_05.travelled_ourDLC, 10*5)*5/100;
dlc_travelled_jc20181219_2_05.travelled_ourDLC = travelled_ourDLC_5fps;

load data\travelled\pix_travelled_jc20181211_1_05.mat
pix_travelled_jc20181211_1_05 = pixels_travelled;
travelled_pixels = medfilt1(pix_travelled_jc20181211_1_05.travelled_pixels, 10*5);
pix_travelled_jc20181211_1_05.travelled_pixels = travelled_pixels;

load data\travelled\pix_travelled_jc20181218_1_05.mat
pix_travelled_jc20181218_1_05 = pixels_travelled;
travelled_pixels = medfilt1(pix_travelled_jc20181218_1_05.travelled_pixels, 10*5);
pix_travelled_jc20181218_1_05.travelled_pixels = travelled_pixels;

load data\travelled\pix_travelled_jc20181219_2_05.mat
pix_travelled_jc20181219_2_05 = pixels_travelled;   
travelled_pixels = medfilt1(pix_travelled_jc20181219_2_05.travelled_pixels, 10*5);
pix_travelled_jc20181219_2_05.travelled_pixels = travelled_pixels;

load data\jc20181211_1\jc20181211_1-seizureOnsetsDatenum.mat
jc20181211_1_seizures = seizureOnsetsDatenum;
load data\jc20181218_1\jc20181218_1-seizureOnsetsDatenum.mat
jc20181218_1_seizures = seizureOnsetsDatenum;
load data\jc20181219_2\jc20181219_2-seizureOnsetsDatenum.mat
jc20181219_2_seizures = seizureOnsetsDatenum;

%% Date of birth each mouse
% subtract day of birth 'jc20181211_1'	
date_of_birth = '20181030';
Y = str2double(date_of_birth(1:4));
M = str2double(date_of_birth(5:6));
D = str2double(date_of_birth(7:8));
jc20181211_1_dob = datenum(datetime(Y,M,D));

% subtract day of birth 'jc20181218_1'	
date_of_birth = '20181104';
Y = str2double(date_of_birth(1:4));
M = str2double(date_of_birth(5:6));
D = str2double(date_of_birth(7:8));
jc20181218_1_dob = datenum(datetime(Y,M,D));

% subtract day of birth 'jc20181219_2'	
date_of_birth = '20181101';
Y = str2double(date_of_birth(1:4));
M = str2double(date_of_birth(5:6));
D = str2double(date_of_birth(7:8));
jc20181219_2_dob = datenum(datetime(Y,M,D));

jc20181211_cut = (dlc_travelled_jc20181211_1_05.t_dlc - jc20181211_1_dob) <= 70;
jc20181211_cut_pix = (pix_travelled_jc20181211_1_05.t_pix - jc20181211_1_dob) <= 70;

dlc_travelled_jc20181211_1_05.travelled_ourDLC = dlc_travelled_jc20181211_1_05.travelled_ourDLC(jc20181211_cut);
dlc_travelled_jc20181211_1_05.t_dlc = dlc_travelled_jc20181211_1_05.t_dlc(jc20181211_cut);

pix_travelled_jc20181211_1_05.travelled_pixels = pix_travelled_jc20181211_1_05.travelled_pixels(jc20181211_cut_pix);
pix_travelled_jc20181211_1_05.t_pix = pix_travelled_jc20181211_1_05.t_pix(jc20181211_cut_pix);

%%

fs = 5; % fps

n_jc20181211_dlc = length(dlc_travelled_jc20181211_1_05.travelled_ourDLC);
f_jc20181211_dlc = (0:n_jc20181211_dlc-1)*(fs/n_jc20181211_dlc);
n_jc20181218_dlc = length(dlc_travelled_jc20181218_1_05.travelled_ourDLC);
f_jc20181218_dlc = (0:n_jc20181218_dlc-1)*(fs/n_jc20181218_dlc);
n_jc20181219_dlc = length(dlc_travelled_jc20181219_2_05.travelled_ourDLC);
f_jc20181219_dlc = (0:n_jc20181219_dlc-1)*(fs/n_jc20181219_dlc);

n_jc20181211_pix = length(pix_travelled_jc20181211_1_05.travelled_pixels);
f_jc20181211_pix = (0:n_jc20181211_pix-1)*(fs/n_jc20181211_pix);
n_jc20181218_pix = length(pix_travelled_jc20181218_1_05.travelled_pixels);
f_jc20181218_pix = (0:n_jc20181218_pix-1)*(fs/n_jc20181218_pix);
n_jc20181219_pix = length(pix_travelled_jc20181219_2_05.travelled_pixels);
f_jc20181219_pix = (0:n_jc20181219_pix-1)*(fs/n_jc20181219_pix);

%% FFT

x1 = dlc_travelled_jc20181211_1_05.travelled_ourDLC - mean(dlc_travelled_jc20181211_1_05.travelled_ourDLC);
y_jc20181211_dlc = fft(x1);
power_jc20181211_dlc = abs(y_jc20181211_dlc).^2/n_jc20181211_dlc;
x2 = dlc_travelled_jc20181218_1_05.travelled_ourDLC - mean(dlc_travelled_jc20181218_1_05.travelled_ourDLC);
y_jc20181218_dlc = fft(x2);
power_jc20181218_dlc = abs(y_jc20181218_dlc).^2/n_jc20181218_dlc;
x3 = dlc_travelled_jc20181219_2_05.travelled_ourDLC - mean(dlc_travelled_jc20181219_2_05.travelled_ourDLC);
y_jc20181219_dlc = fft(x3);
power_jc20181219_dlc = abs(y_jc20181219_dlc).^2/n_jc20181219_dlc;

x4 = pix_travelled_jc20181211_1_05.travelled_pixels - mean(pix_travelled_jc20181211_1_05.travelled_pixels);
y_jc20181211_pix = fft(x4);
power_jc20181211_pix = abs(y_jc20181211_pix).^2/n_jc20181211_pix;
x5 = pix_travelled_jc20181218_1_05.travelled_pixels - mean(pix_travelled_jc20181218_1_05.travelled_pixels);
y_jc20181218_pix = fft(x5);
power_jc20181218_pix = abs(y_jc20181218_pix).^2/n_jc20181218_pix;
x6 = pix_travelled_jc20181219_2_05.travelled_pixels - mean(pix_travelled_jc20181219_2_05.travelled_pixels);
y_jc20181219_pix = fft(x6);
power_jc20181219_pix = abs(y_jc20181219_pix).^2/n_jc20181219_pix;

%% visualise
close all;

% figure(10);
% 
% subplot(311);
% pwelch(x1);
% 
% subplot(312);
% pwelch(x2);
% 
% subplot(313);
% pwelch(x3);
% 
% figure(11);
% 
% subplot(311);
% periodogram(x1);
% 
% subplot(312);
% periodogram(x2);
% 
% subplot(313);
% periodogram(x3);

[max1, i1] = max(power_jc20181211_dlc);
[max2, i2] = max(power_jc20181218_dlc);
[max3, i3] = max(power_jc20181219_dlc);

maxf1 = f_jc20181211_dlc(i1);
maxf2 = f_jc20181218_dlc(i2);
maxf3 = f_jc20181219_dlc(i3);

maxt1 = 1/maxf1 /60/60/24;
maxt2 = 1/maxf2 /60/60/24;
maxt3 = 1/maxf3 /60/60/24;

max_days = 10;
Xlim = [0 max_days];
Xlabel = 'Time (day)';
Xticks = 0:max_days;
Ylabel = 'Power ((m/s)^2/Hz)';

figure(10);
plot(1./f_jc20181211_dlc/60/60/24, power_jc20181211_dlc, 'k');
title('DLC jc20181211\_1');
xlim(Xlim);
xlabel(Xlabel);
xticks(Xticks);
ylabel(Ylabel);
fontsize(gcf,scale=3);

figure(11);
plot(1./f_jc20181218_dlc/60/60/24, power_jc20181218_dlc, 'k');
title('DLC jc20181218\_1');
xlim(Xlim);
xlabel(Xlabel);
xticks(Xticks);
ylabel(Ylabel);
fontsize(gcf,scale=3);

figure(12);
plot(1./f_jc20181219_dlc/60/60/24, power_jc20181219_dlc, 'k');
title('DLC jc20181219\_2');
xlim(Xlim);
xlabel(Xlabel);
xticks(Xticks);
ylabel(Ylabel);
fontsize(gcf,scale=3);

Ylabel = 'Power (-)';
figure(13);
plot(1./f_jc20181211_pix/60/60/24, power_jc20181211_pix, 'k');
title('jc20181211\_1');
xlim(Xlim);
xlabel(Xlabel);
xticks(Xticks);
ylabel(Ylabel);
fontsize(gcf,scale=3);

figure(14);
plot(1./f_jc20181218_pix/60/60/24, power_jc20181218_pix, 'k');
title('jc20181218\_1');
xlim(Xlim);
xlabel(Xlabel);
xticks(Xticks);
ylabel(Ylabel);
fontsize(gcf,scale=3);

figure(15);
plot(1./f_jc20181219_pix/60/60/24, power_jc20181219_pix, 'k');
title('jc20181219\_2');
xlim(Xlim);
xlabel(Xlabel);
xticks(Xticks);
ylabel(Ylabel);
fontsize(gcf,scale=3);
