clear all
close all

load('transient_oscillations.mat')

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

final_=final(1:n/2)
t_ = t(1:n/2)
n_=n/2

figure(1)
subplot(3,1,1)
plot(t_,final_)
xlabel('time')
ylabel('data')

subplot(3,1,2)
% compute FT
num_windows=3;
num_window_points=n_/num_windows
duration_window=duration/num_windows
fs=1/dt;
overlap=num_window_points-1
%window = rectwin(num_window_points);
window = hamming(num_window_points);
[S,F,T,P]=spectrogram(final_,window,overlap,num_window_points,fs);
%size(P)
%T
%Pt=P'
pcolor(T,F,P);
%plot(10*log(Pt(1,:)))
Pmax=max(max(P))
Pmin=min(min(P))
Tmax=max(max(T))
Tmin=min(min(T))
view(0,90);
caxis([0,Pmax]);shading flat;
axis([Tmin,Tmax,0,2.0]);
colorbar();
ylabel('frequency')
xlabel('time [s]')


% 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,3)
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,'scales',scales_opts);
    PP(:,i)=abs(cwtstruct.cfs).^2;
end
logPP=10*log10(PP);
minPP=min(squeeze(min(PP)))
maxPP=max(squeeze(max(PP)))
minlogPP=min(squeeze(min(logPP)));
maxlogPP=max(squeeze(max(logPP)));
%pcolor(squeeze(time),freq,logPP');
pcolor(squeeze(t_),freq,PP');
% %imagesc(squeeze(data_1(:,1)),freq,10*log10(PP));
view(0,90);
%caxis([minlogPP,maxlogPP]);shading flat;
caxis([minPP,maxPP]);shading flat;
%caxis([20,maxlogPP]);shading flat;
%caxis([1.0e+02,2.0e+06]);shading flat;
axis([T0,T1,fmin,fmax-df]);
colorbar();
ylabel('frequency')
xlabel('time [s]')


