% Script written by Tanya Leise, Department of Mathematics, Amherst College
% Calls scripts wavetrans.m and ridgewalk.m from jlab written by J. M. Lilly
% (http://www.jmlilly.net/jmlsoft.html)
%
% Also uses custom scripts Actogram.m, CalcScaleForCWT.m, and CWTPlotLabels.m
%
% This script plots the analytic wavelet transform of a mouse wheel-running record.
% Also calculates and plots the wavelet ridge curves.


clear

longestperiod=30;
shortestperiod=20;
gamma=3;beta=7;  % always use gamma=3, decrease beta to increase time resolution (at expense of frequency resolution)

filename = 'OUTSYNC 170 Wheel-TimeSeries.csv';
M=dlmread(filename,',',1,0);
t=M(:,1);
Nsteps=numel(t);
wheel=M(:,2);
light=M(:,3);
clear M

NstepsPerHr=4;TimeStep=1/NstepsPerHr;  % wheel revolutions are in 15 minute bins
T=t(Nsteps)-t(1);

Actogram(t,wheel,light,'new',24)
title(filename);

% calculate wavelet transform and ridge
nvoice=512;
[fs,tau,qscaleArray] = CalcScaleForCWT(shortestperiod,longestperiod,T,NstepsPerHr,nvoice);

cwt=wavetrans(wheel,{1,gamma,beta,fs,'bandpass'},'periodic');  % cwt is a matrix with complex numbers (jlab)

[wr,ir,jr,fr]=ridgewalk(TimeStep,cwt,fs);
keep=find(~isnan(ir) & 2*pi./fr>18 & 2*pi./fr<28);
wr=wr(keep);jr=jr(keep);fr=fr(keep);ir=ir(keep);
tridge=t(ir)/24;
ridgetau=2*pi./fr; % period estimate along ridge
ridgeampl=abs(wr); % amplitude estimate along ridge

%  wr     Wavelet transfrom value along the ridge
%  ir     Ridge indices into rows of cwt (time) 
%  jr     Ridge indices into columns of cwt (scale) 
%  fr     Instantaneous frequency along the ridge

% plot wavelet coefficients (scalogram)
figure('OuterPosition',[100 500 800 400]);
imagesc(abs(cwt)');colorbar; % plots the wavelet transform coefficients, which give the amplitude
hold on;scatter(ir,jr,4,'k','filled')  % plots the ridges

tauticks=20:2:28;
tticks=(24*4):(24*4):t(end);
CWTPlotLabels(nvoice,qscaleArray,'integer',tau,tauticks,t,tticks)
set(gca,'FontSize',20);
xlabel('Time in days')
ylabel('Period (h)');
title('AWT of wheel-running activity')

% plot AWT phase estimate at each time of lights on
dawn=[3.2472:25.2472 25.9139:31.9139 32.5806:62.5806];
dawnind=round((dawn*24-t(1))*NstepsPerHr+1);
dawnir=interp1(ir,1:numel(ir),dawnind);
arg=angle(cwt(ir+Nsteps*(jr-1)));
phase=mod(arg(dawnir)*12/pi,24);

figure('OuterPosition',[50 500 600 450]);
plot(dawn,phase,'k','LineWidth',2)
set(gca,'FontSize',20,'XTick',4:8:60);axis([0 64 0 10]);
xlabel('Time in days');ylabel('Phase (h)')

% plot AWT-estimated period and amplitude over time
figure('OuterPosition',[0 500 600 450]);
yyaxis left
plot(tridge,ridgetau,'b','LineWidth',2)
xlim([42 61.75])

yyaxis right
plot(tridge,ridgeampl,'r-.','LineWidth',2)
xlabel('Time in days')
set(gca,'FontSize',20)
legend('Period','Amplitude','Location','SouthWest')

ax=gca;
ax.YAxis(1).Limits=[22 25];
ax.YAxis(1).TickValues=22:25;
ax.YAxis(2).Limits=[340 510];
ax.YAxis(2).TickValues=350:50:500;
