% 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 CalcScaleForCWT.m and CWTPlotLabels.m
%
% This script plots three different signals and their analytic wavelet
% transforms to illustrate how it works on signals with varying frequency or amplitude. 

clear all

NstepsPerHr=6;
TimeStep=1/NstepsPerHr;

period1=24;  % period in hours of main oscillation
period2=24*7;% period in hours governing the varying frequency or amplitude

Nsteps=24*14*NstepsPerHr+1;
t=(0:Nsteps-1)'/NstepsPerHr;
T=t(Nsteps)-t(1);

longestperiod=42;shortestperiod=14;  % range of periods to use in AWT

gamma=3;beta=7; % parameter values for Morse wavelet function

nvoice=128;
[fs,tau,qscaleArray] = CalcScaleForCWT(shortestperiod,longestperiod,T,NstepsPerHr,nvoice);
taubnds=round(interp1(tau,1:numel(tau),[14 37]));  

% Basic sinusoid
figure;
x=sin(2*pi*t/period1);  
plot(t,x,'k','LineWidth',2);
set(gca,'FontSize',20,'PlotBoxAspectRatio',[1 0.6 1],'Position',[.13 .11 .775 .76]);
axis([0 T -1.1 1.1]);
days=0:48:t(Nsteps);
C=cell(1,numel(days));
for k=1:numel(days)
    C(1,k)={num2str(2*(k-1))};
end
set(gca,'XTickLabel',C,'XTick',days);
xlabel('Time in days')
title('Sinusoid')

cwt=wavetrans(x,{1,gamma,beta,fs,'bandpass'},'periodic');
[ir,jr,~,~]=ridgewalk(cwt,fs);

figure;
imagesc(abs(cwt)');colorbar;caxis([0 1])
hold on;scatter(ir,jr,6,'k','filled')
CWTPlotLabels(t(1),t(Nsteps),NstepsPerHr,nvoice,qscaleArray,'integer',tau,shortestperiod,longestperiod,0,0,2);
set(gca,'FontSize',20,'PlotBoxAspectRatio',[1 0.6 1]);
xlabel('Time in days')
ylim(taubnds);ylabel('Period (h)');xlabel('');
title('AWT of sinusoid')

% Varying frequency
df=1;
x=sin(2*pi*t/period1+df*sin(2*pi*t/period2)); 
figure;
plot(t,x,'k','LineWidth',2);
set(gca,'FontSize',20,'PlotBoxAspectRatio',[1 0.6 1],'Position',[.13 .11 .775 .76]);
axis([0 T -1.1 1.1]);
days=0:48:t(Nsteps);
C=cell(1,numel(days));
for k=1:numel(days)
    C(1,k)={num2str(2*(k-1))};
end
set(gca,'XTickLabel',C,'XTick',days);
xlabel('Time in days')
title('Signal with frequency variation')

cwt=wavetrans(x,{1,gamma,beta,fs,'bandpass'},'periodic');
[ir,jr,~,~]=ridgewalk(cwt,fs);

figure;
imagesc(abs(cwt)');colorbar;caxis([0 1])
hold on;scatter(ir,jr,6,'k','filled')
CWTPlotLabels(t(1),t(Nsteps),NstepsPerHr,nvoice,qscaleArray,'integer',tau,shortestperiod,longestperiod,0,0,2);
set(gca,'FontSize',20,'PlotBoxAspectRatio',[1 0.6 1]);
xlabel('Time in days')
ylim(taubnds);ylabel('Period (h)');xlabel('');
title('AWT with varying frequency ')

% Varying amplitude
x=sin(2*pi*t/period1).*(.8+.2*sin(2*pi*t/period2));
figure;
plot(t,x,'k','LineWidth',2);
set(gca,'FontSize',20,'PlotBoxAspectRatio',[1 0.6 1],'Position',[.13 .11 .775 .76]);
axis([0 T -1.1 1.1]);
days=0:48:t(Nsteps);
C=cell(1,numel(days));
for k=1:numel(days)
    C(1,k)={num2str(2*(k-1))};
end
set(gca,'XTickLabel',C,'XTick',days);
xlabel('Time in days')
title('Signal with varying amplitude ')

cwt=wavetrans(x,{1,gamma,beta,fs,'bandpass'},'periodic');
[ir,jr,~,~]=ridgewalk(cwt,fs);

figure;
imagesc(abs(cwt)');colorbar;caxis([0 1])
hold on;scatter(ir,jr,6,'k','filled')
CWTPlotLabels(t(1),t(Nsteps),NstepsPerHr,nvoice,qscaleArray,'integer',tau,shortestperiod,longestperiod,0,0,2);
set(gca,'FontSize',20,'PlotBoxAspectRatio',[1 0.6 1]);

xlabel('Time in days')
ylim(taubnds);ylabel('Period (h)');xlabel('');
title('AWT with amplitude variation')