
%% load the data
data=load('yearly_outside_temperatures.dat');
t=data(:,1);   % time in seconds (unix date format)
d=t/24/60/60;  % time conversion to day number convention
date_offset=datenum([1970 1 1 0 0 t(1)]); % unix epoch to Matlab
d=d-min(d);    % removing large unnecessary offset
d=d+date_offset; % converting to Matlab date format
to=data(:,2);  % outside temperature
Ndata=length(to);  % how many points

figure(1);
plot(d, to, '.');
ylabel('Temp (^oC)');
title('Raw data');
datetick('x','yyyy/mmm'); % x to date

%% Have a look at the temperature data array length 
data_length=length(to)

%% Average spacing in seconds between consequent points excluding longer then hour voids
dt=diff(t); 
data_points_spacing=mean(dt(dt <= 3600))

%% let's reduce the data size for now with running average
% the size of yearly variation  is about 10'000 points
% lets reduce the data size, 
% besides we do not need the every minute data most likely it is due to clouds and rain
Nf=60*6; % designed to be one data point per 6 hours
dr=reduce_with_average(d,Nf);
tor=reduce_with_average(to,Nf);
figure(1)
plot(dr, tor);
ylabel('Temp (^oC)');
datetick('x','yyyy/mmm'); % x to date
title('Reduced data set');

%% There are a few data gaps. One especially big is in 2014.
% we will do BIG NO-NO and fill the gap with the interpolated data
tstart=dr(1); tend=dr(end);

% setting sample spacing
dt=1/4; % recall that matlab date has number of day before decimal
t=tstart:1/4:tend; % additional benefit of fixing a distance between samples
Temp_Interpolated=interp1(dr,tor, t);

figure(1)
plot(t, Temp_Interpolated);
ylabel('Temp (^oC)');
datetick('x','yyyy/mmm'); % x to date
title('Interpolated data set');

%% now let' see the Fourier spectrum of the weather data
Y=fft(Temp_Interpolated-mean(Temp_Interpolated) ); % we get rid of mean temperature
N=length(Temp_Interpolated);
Y=Y/N; % normalizing to amplitudes of sin/cos since fft is missing this factor 1/N

freq=fourier_frequencies(1/dt, N);

figure(2);
plot(freq, abs(Y), '.');
xlabel('Frequency (1/day)');
title('Spectrum of the temperature data set');

%% since we plot only |Y| there is no need to show negative frequency 
% since they have the same amplitude
Y=Y(freq>0);
freq=freq(freq>0);

figure(2);
plot(freq, abs(Y), '-');
xlabel('Frequency (1/day)');
title('Spectrum of the temperature data set');

%% let's convert frequency to period
period=1./freq;

figure(3);
semilogx(period, abs(Y), '-');
xlabel('Oscillation period  (day)');
title('Amplitude of oscillation of the temperature data set');


