clear all
close all

%load('transient_oscillations.mat')

% generate data example
duration = 20;
duration_activity=5;
activity_start=(duration-duration_activity)/2;
dt=0.02;
n = duration / dt
fmax=1/dt;

t = linspace(0,duration-dt,n);
final=zeros(n,1);
n_activity_start = activity_start/dt;
n_activity=duration_activity/dt;
for i = 0:n_activity-1
    final(n_activity_start+i)=sin(2*pi*1.0*dt*i) ;   
end

figure(1)
subplot(3,1,1)
plot(t,final);
xlabel('time')
ylabel('data')

subplot(3,1,2)
% compute FT
num_windows=1;
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=4;
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=3;
pxx4 = pmtm(final,nw,n);
nw=5;
pxx5 = pmtm(final,nw,n);

% %plot(f1(1:num_maxf1),pxx1(1:num_maxf1),f2(1:num_maxf2),pxx2(1:num_maxf2),'r')
% %plot(f1(1:num_maxf1),pxx1(1:num_maxf1))
% xlabel('frequency')
% ylabel('power spectral density [dB]')


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

subplot(3,1,3)
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]')
