clear all
close all

load('transient_oscillations.mat')

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

fs=1/dt
order = 4

lowFreq=0.9
hiFreq=1.1
[b,a] = butter(order, [lowFreq hiFreq]/(fs/2), 'bandpass');
 y1 = filtfilt(b,a,final)

lowFreq=0.5
hiFreq=0.7
[b,a] = butter(order, [lowFreq hiFreq]/(fs/2), 'bandpass');
 y2 = filtfilt(b,a,final)
 
lowFreq=0.3
hiFreq=0.5
[b,a] = butter(order, [lowFreq hiFreq]/(fs/2), 'bandpass');
 y3 = filtfilt(b,a,final)

 figure(1)
subplot(4,1,1)
plot(t,final)
xlabel('time')
ylabel('signal')
subplot(4,1,2)
plot(t,y1)
xlabel('time')
ylabel('band 1Hz')
subplot(4,1,3)
plot(t,y2)
xlabel('time')
ylabel('band 0.6Hz')
subplot(4,1,4)
plot(t,y3)
xlabel('time')
ylabel('band 0.4Hz')

[r1,lags1] = xcorr(y1,'coeff'); 
[r2,lags2] = xcorr(y2,'coeff'); 
[r3,lags3] = xcorr(y3,'coeff'); 

figure(2)
subplot(3,1,1)
plot(lags1*dt,r1)
xlabel('time')
ylabel('band 1Hz')
subplot(3,1,2)
plot(lags2*dt,r2)
xlabel('time')
ylabel('band 0.6Hz')
subplot(3,1,3)
plot(lags3*dt,r3)
xlabel('time')
ylabel('band 0.4Hz')

meanFreq=1.0

dist=0.2
lowFreq=meanFreq-dist
hiFreq=meanFreq+dist
[b,a] = butter(order, [lowFreq hiFreq]/(fs/2), 'bandpass');
 y4 = filtfilt(b,a,final)

dist=0.3
lowFreq=meanFreq-dist
hiFreq=meanFreq+dist
[b,a] = butter(order, [lowFreq hiFreq]/(fs/2), 'bandpass');
 y5 = filtfilt(b,a,final)
 
figure(3)
subplot(3,1,1)
plot(t,y1)
xlabel('time')
ylabel('bandwidth 0.2Hz')
subplot(3,1,2)
plot(t,y4)
xlabel('time')
ylabel('bandwidth 0.4Hz')
subplot(3,1,3)
plot(t,y5)
xlabel('time')
ylabel('bandwidth 0.6Hz')

[r4,lags4] = xcorr(y4,'coeff'); 
[r5,lags5] = xcorr(y5,'coeff'); 

figure(5)
subplot(3,1,1)
plot(lags1*dt,r1)
xlabel('time')
ylabel('bandwidth 0.2Hz')
subplot(3,1,2)
plot(lags4*dt,r4)
xlabel('time')
ylabel('bandwidth 0.4Hz')
subplot(3,1,3)
plot(lags5*dt,r5)
xlabel('time')
ylabel('bandwidth 0.6Hz')