% 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

NstepsPerHr=6; % # samples per hour
TimeStep=1/NstepsPerHr; % interval between time points

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; % total # of time points
t=(0:Nsteps-1)'/NstepsPerHr; % time vector
T=t(Nsteps)-t(1); % total duration

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

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

nvoice=64; % resolution for AWT coefficients
[fs,tau,qscaleArray] = CalcScaleForCWT(shortestperiod,longestperiod,T,NstepsPerHr,nvoice);
taubnds=round(interp1(tau,1:numel(tau),[14 37]));  

% Basic sinusoid
x=cos(2*pi*t/period1);  

figure('Position',[1 472 900 350]);
plot(t/24,x,'k','LineWidth',2);
set(gca,'FontSize',20,'XTick',0:T/24);
axis([0 T/24 -1.1 1.1]);
xlabel('Time in days')
title('Simple sinusoid')

cwt=wavetrans(x,{1,gamma,beta,fs,'bandpass'},'mirror');
[wr,ir,jr,fr]=ridgewalk(TimeStep,cwt,fs);
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

figure('Position',[1 472 900 350]);
imagesc(abs(cwt)');colorbar;caxis([0 1])
hold on;scatter(ir,jr,6,'k','filled')

tauticks=16:4:40;
tticks=0:24:T;
CWTPlotLabels(nvoice,qscaleArray,'integer',tau,tauticks,t,tticks)
set(gca,'FontSize',20);
xlabel('Time in days')
ylim(taubnds);ylabel('Period (h)');
title('AWT of simple sinusoid')

figure('Position',[1 472 900 350]);
indcirc=find(abs(ridgetau-24)<4);
plot(t(ir(indcirc))/24,ridgeampl(indcirc))
set(gca,'FontSize',20,'ylim',[0 1.5],'xtick',0:T/24);
xlabel('Time in days')
ylabel('Amplitude');
title('Simple sinusoid example')

% Varying frequency
df=1;
x=cos(2*pi*t/period1+df*sin(2*pi*t/period2)); 

figure('Position',[1 472 900 350]);
plot(t/24,x,'k','LineWidth',2);
set(gca,'FontSize',20,'XTick',0:T/24);
axis([0 T/24 -1.1 1.1]);
xlabel('Time in days')
title('Signal with frequency varying over time')

cwt=wavetrans(x,{1,gamma,beta,fs,'bandpass'},'mirror');
[wr,ir,jr,fr]=ridgewalk(TimeStep,cwt,fs,1,2); % extra 1,2 suppresses short artifactual ridges
ridgetau=2*pi./fr; % period estimate along ridge
ridgeampl=abs(wr); % amplitude estimate along ridge

figure('Position',[1 472 900 350]);
imagesc(abs(cwt)');colorbar;caxis([0 1])
hold on;scatter(ir,jr,6,'k','filled')

tauticks=16:4:40;
tticks=0:24:T;
CWTPlotLabels(nvoice,qscaleArray,'integer',tau,tauticks,t,tticks)
set(gca,'FontSize',20);
xlabel('Time in days')
ylim(taubnds);ylabel('Period (h)');
title('AWT for varying frequency example')

figure('Position',[1 472 900 350]);
indcirc=find(abs(ridgetau-24)<4);
plot(t(ir(indcirc))/24,ridgeampl(indcirc))
set(gca,'FontSize',20,'ylim',[0 1.5],'xtick',0:T/24);
xlabel('Time in days')
ylabel('Amplitude');
title('Varying frequency example')

% Varying amplitude
x=cos(2*pi*t/period1).*(.8+.2*cos(2*pi*t/period2));

figure('Position',[1 472 900 350]);
plot(t/24,x,'k','LineWidth',2);
set(gca,'FontSize',20,'XTick',0:T/24);
axis([0 T/24 -1.1 1.1]);
xlabel('Time in days')
title('Signal with varying amplitude')

cwt=wavetrans(x,{1,gamma,beta,fs,'bandpass'},'mirror');
[wr,ir,jr,fr]=ridgewalk(TimeStep,cwt,fs);
ridgetau=2*pi./fr; % period estimate along ridge
ridgeampl=abs(wr); % amplitude estimate along ridge

figure('Position',[1 472 900 350]);
imagesc(abs(cwt)');colorbar;caxis([0 1])
hold on;scatter(ir,jr,6,'k','filled')

tauticks=16:4:40;
tticks=0:24:T;
CWTPlotLabels(nvoice,qscaleArray,'integer',tau,tauticks,t,tticks)
set(gca,'FontSize',20);
xlabel('Time in days')
ylim(taubnds);ylabel('Period (h)');
title('AWT for varying amplitude example')

figure('Position',[1 472 900 350]);
indcirc=find(abs(ridgetau-24)<4);
plot(t(ir(indcirc))/24,ridgeampl(indcirc))
set(gca,'FontSize',20,'ylim',[0 1.5],'xtick',0:T/24);
xlabel('Time in days')
ylabel('Amplitude');
title('Varying amplitude example')
