clear all
close all

load('transient_oscillations.mat')

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

figure(1)
% subplot(2,1,1)
% plot(t,final)
% xlabel('time')
% ylabel('data')

subplot(2,1,1)
% compute FT
% num_windows=5;
% num_window_points=n/num_windows;
% num_overlap = 0.0 * num_window_points;
% duration_window=duration/num_windows;
% df=1/duration_window;
% num_maxf1=2.0/df;
% [pxx1,f1] = pwelch(final,num_window_points,num_overlap,num_window_points,fmax);

num_windows=5;
num_window_points=n/num_windows;
num_overlap = 0.5 * num_window_points;
duration_window=duration/num_windows;
df=1/duration_window;
[pxx2,f2] = pwelch(final,num_window_points,num_overlap,num_window_points,fmax);
%plot(f,10*log10(pxx))
num_maxf2=2.0/df;

df=1/duration;
num_maxf3=2.0/df;
freq3 = linspace(0,fmax/2,n/2);
nw=2;
pxx3 = pmtm(final,nw,n);
nw=5;
pxx4 = pmtm(final,nw,n);
nw=10;
pxx5 = pmtm(final,nw,n);

%plot(f1(1:num_maxf1),pxx1(1:num_maxf1),f2(1:num_maxf2),pxx2(1:num_maxf2),'r')
plot(f2(1:num_maxf2),10*log(pxx2(1:num_maxf2)),'r',freq3(1:num_maxf3),10*log(pxx3(1:num_maxf3)),'g')
%plot(f1(1:num_maxf1),pxx1(1:num_maxf1))
xlabel('frequency')
ylabel('power spectral density [dB]')

subplot(2,1,2)
plot(freq3(1:num_maxf3),10*log(pxx3(1:num_maxf3)),'g',freq3(1:num_maxf3),10*log(pxx4(1:num_maxf3)),'m',freq3(1:num_maxf3),10*log(pxx5(1:num_maxf3)),'b')
%plot(f1(1:num_maxf1),pxx1(1:num_maxf1))
xlabel('frequency')
ylabel('power spectral density [dB]')
