clear all
close all

load('transient_oscillations.mat')

% generate data example
duration = max(t)
n = duration / dt+1


% wavelet power spectrum
t0=0;
T0=0.0;
t1=n;% number of points
T1=n*dt; % in seconds
fmin=0.0;
fmax=2.0;
num_scales=200;
df=(fmax-fmin)/num_scales;
MorletFourierFactor = 4.0*pi/(6.0+sqrt(2.0+6.0^2));
nb=1

subplot(3,1,1)
sig = struct('val',final,'period',dt);
sum_alpha=0;
alpha_counter=0;
for i = 1:num_scales
    scale = i-1;
    freq(i)=fmin+scale*df;
    scale_low = 1.0/(MorletFourierFactor*freq(i));
    ds = 1.0/(MorletFourierFactor*(fmin+scale*df))-1.0/(MorletFourierFactor*(fmin+(scale+1)*df));
    scales_opts = struct('s0',scale_low,'ds',ds,'nb',nb,'type','lin');
    cwtstruct = cwtft(sig,'wavelet','morl','scales',scales_opts);
    PP(:,i)=abs(cwtstruct.cfs).^2;
end
minPP=min(squeeze(min(PP)))
maxPP=max(squeeze(max(PP)))
pcolor(squeeze(t),freq,PP');
view(0,90);
caxis([minPP,maxPP]);shading flat;
axis([T0,T1,fmin,fmax-df]);
colorbar();
ylabel('frequency')
xlabel('time [s]')

subplot(3,1,2)
sum_alpha=0;
alpha_counter=0;
for i = 1:num_scales
    scale = i-1;
    freq(i)=fmin+scale*df;
    scale_low = 1.0/(MorletFourierFactor*freq(i));
    ds = 1.0/(MorletFourierFactor*(fmin+scale*df))-1.0/(MorletFourierFactor*(fmin+(scale+1)*df));
    scales_opts = struct('s0',scale_low,'ds',ds,'nb',nb,'type','lin');
    cwtstruct = cwtft(sig,'wavelet','morlex','scales',scales_opts);
    PP(:,i)=abs(cwtstruct.cfs).^2;
end
minPP=min(squeeze(min(PP)))
maxPP=max(squeeze(max(PP)))
pcolor(squeeze(t),freq,PP');
view(0,90);
caxis([minPP,maxPP]);shading flat;
axis([T0,T1,fmin,fmax-df]);
colorbar();
ylabel('frequency')
xlabel('time [s]')

subplot(3,1,3)
sum_alpha=0;
alpha_counter=0;
for i = 1:num_scales
    scale = i-1;
    freq(i)=fmin+scale*df;
    scale_low = 1.0/(MorletFourierFactor*freq(i));
    ds = 1.0/(MorletFourierFactor*(fmin+scale*df))-1.0/(MorletFourierFactor*(fmin+(scale+1)*df));
    scales_opts = struct('s0',scale_low,'ds',ds,'nb',nb,'type','lin');
    cwtstruct = cwtft(sig,'wavelet','morl0','scales',scales_opts);
    PP(:,i)=abs(cwtstruct.cfs).^2;
end
minPP=min(squeeze(min(PP)))
maxPP=max(squeeze(max(PP)))
pcolor(squeeze(t),freq,PP');
view(0,90);
caxis([minPP,maxPP]);shading flat;
axis([T0,T1,fmin,fmax-df]);
colorbar();
ylabel('frequency')
xlabel('time [s]')