%% Frederik J Simons - Lab04 - fjsimonsl04
%% Read the data
% The usual: webread. It's quite a big file that we made for you, so we're loading it in a Section all by itself.

w=webread('http://geoweb.princeton.edu/people/simons/weather_data.csv');

% Make some strings up on the top here, save us the typing later on.

ex1='calendar day';
ex2='hours';
wi1=sprintf('outside temperature (%sC)',char(176));
wot='Hargraves Hall';
ex3='frequency in hours^{-1}';
wi2='variance of the temperature per frequency';
wot2=sprintf('periodogram of the temperature at %s',wot);

%
%
%% Make a simple plot using calendar time
% There is a bit of overhead here, related to how our weather station stores time: namely, in Unix timestamps, and that means: in seconds since 1/1/1970 00:00:00 GMT. But we will use MATLABs datenum, datetick and datestr to deal with these types of calendar data. Be prepared to go back and forth! Also, of course, you need to just "know" that column 1 contains the time and column 5 the outdoor temperature data collected on the roof of Hargraves Hall in the Fall of 2015.

% Use the explicit structure addressing that gets created when a CSV table has no column headers
unixtime=w.Var1;
temperature=w.Var5;
% Convert the date format into MATLABs own date format, which is in days, not seconds
matlabtime=unixtime/60/60/24+datenum(1970,1,1,0,0,0);
% Make a simple plot
subplot(211)
plot(matlabtime,temperature); axis tight
% Use this trick to set the ticks to calendar dates
datetick('x',6)
% Label, title, try to be as informative as possible
xlabel(ex1); ylabel(wi1)
title(sprintf('%s between %s and %s',...
    wot,datestr(matlabtime(end)),datestr(matlabtime(1)))) 
grid on 

%
%
%% Make a simple plot using hourly time
% Our objective will be to determine the periodicity that we suspect is present in this time series, so we will replot the data with a horizontal axis of hours this time. Do not call the variable hours, which is a MATLAB function itself.

% This line is the only thing that really happens in this section
mlhours=[matlabtime-matlabtime(end)]*24;
% I'm preparing a "tick" mark multiplier - see below
hmult=5;
% All the rest is merely replotting as before
subplot(211)
plot(mlhours,temperature); axis tight
set(gca,'xtick',0:24*hmult:mlhours(1))
xlabel(sprintf('%s between %s and %s',ex2,...
    datestr(matlabtime(end),'dd mmm yyyy, HH:MM:SS'),...
    datestr(matlabtime(1),'dd mmm yyyy, HH:MM:SS')))
ylabel(wi1); title(wot); grid on 

%
%% Let the data tell us what a fitting sinusoid is
% At the most basic level of our understanding of weather, wouldn't we expect there to be a diurnal cycle? How about, indeed, some kind of a sinusoid above a certain offset, with some kind of an amplitude, a wavelength of, I don't know, 24 hours, and a phase offset that we might need to find by trial and error.
% What would a good starter offset and amplitude be?

% Starter offset guess
m=mean(temperature);
% Starter amplitude guess
A=std(temperature);

% Let's make a sinusoid with these parameters, and also with a wavelength of T hours, and a phase offset of phi, and overlay it on the plot. As a metric for how close we get in fitting the single-sinusoid character of the temperature, play with the period T and with the value phi until you get maximum "agreement" - as measured by the correlation coefficient.

% Make the wavelength and the phase offset a variable also
T=24;
phi=2.3693;
% Attempt to capture some features of the temperature cycle by a single sinusoid
s=m+A*sin(2*pi/T*mlhours+phi);
% Calculate the correlation coefficient between the observed data (temperature) and the modeled data (s)
cst=corr(s(:),temperature(:));
% Overlay on the plot and show the value of the correlation coefficient in the title
try; delete(p); end; hold on
p=plot(mlhours,s,'r');
hold off 
title(sprintf('%s with T %g h and offset %i h %i m for corr %g',...
    wot,T,floor(phi/2/pi*T),round(rem(phi/2/pi*T,floor(phi/2/pi*T))*60),cst))

%
%% In-class exercise: Are you pleased? Sure you cannot do any better?
% What is the best you can do? Try a few periods that all play at the same time!

% Some periods that you think maybe at play
T=[24 1 1];
% Some amplitudes that look like they're working
A=[8 0  0];
% Some phases that you can go with
phi=[2.3 0 0];
% And now, a superposition signal that you, again, compare to the record
s=m+A(1)*sin(2*pi/T(1)*mlhours+phi(1))+...
    A(2)*sin(2*pi/T(2)*mlhours+phi(2))+...
    A(3)*sin(2*pi/T(3)*mlhours+phi(3));
% Calculate the correlation coefficient between the observed data (temperature) and the modeled data (s)
cst=corr(s(:),temperature(:));
% Overlay on the plot and show the value of the correlation coefficient in the title
try; delete(p); end; hold on
p=plot(mlhours,s,'r');
hold off 
% Focus on the report of the correlation coefficient as you rerun the figure
title(sprintf('corr %g',cst))

%
%% The professional way: Fourier analysis
% An automated way of recovering cycle periods and their amplitude is by breaking the signal down into a Fourier series,  and then to plot the magnitude of the (sine and cosine) coefficients on a suitable grid of periods. Here goes, using the built-in functions fft (the Fast Fourier Transform) and hanning (a window that tapers the signal up at the beginning of the series, and down at the end). The result is known as the (windowed) periodogram.

% An automated way to find any and all peaks in these data
nfft=length(mlhours);
physl=mlhours(1);
% Make a frequency axis ahead of time
% Here is the sampling interval
xsint=physl/(length(temperature)-1);
% Here is the set of frequencies that will be "enough" to look at
if rem(nfft,2); selekt=[1:(nfft+1)/2]; else; selekt=[1:nfft/2+1]; end
% Here are the frequencies
fax=(selekt-1)'/xsint/nfft;

% Compute a naive spectral estimate, the windowed periodogram
S=abs(fft(hanning(length(temperature)).*[temperature-mean(temperature)],nfft)).^2;

% Plot the periodogram
clf
semilogx(fax,S(selekt),'linew',2)
grid on 
xlabel(ex3); ylabel(wi2); title(wot2)

%
%
%
%% Quote the periods as an overlay:
% We should be trying to read amplitudes and periods off the Fourier plot, and overlay them on the other one.

% Don't overwrite, but overlay!
hold on
% Two periods that hold particular interest
P=[24 12];
% Plot them as vertical lines
plot(1./[P(:) P(:)]',repmat(ylim,length(P),1)','k')
% Done overlaying!
hold off

%
%% Experiment with a log-log plot
% There is a large dynamic range, so a logarithmic transformation will be handy. Let's redo the plot and its overlay.

% Rather than taking the logs ourselves, we can ask MATLAB to make a log-log plot!
loglog(fax,S(selekt))
% Labels and gridding
xlabel(ex3); ylabel(wi2); title(wot2); grid on
% The overlay, again
hold on; plot(1./[P(:) P(:)]',repmat(ylim,length(P),1)','k'); hold off

%% Fit a regression line to the log-log spectrum
% We see that the spectrum, in log-log-space, drops off almost linerarly. Let's determine what line fits best, using polyfit and polyval as we have done before.

% Compute the regression line
regr=polyfit(log(fax(2:end)),log(S(selekt(2:end))),1);
slope=regr(1);
intrc=regr(2);
% Evaluate the regression line between the end points
y=polyval([slope intrc],log(fax));

% Replot but this time more explicitly
plot(log(fax),log(S(selekt)))
hold on
plot(log(fax),y,'r','linew',2)
xlim([log(fax(2)) log(fax(end))])
hold off
% Labels and gridding - make sure to collect an out handle!
xx=xlabel(ex3); yy=ylabel(wi2); tt=title(wot2); grid on
% The overlay, again, but feed it to logs of the periods!
hold on; plot(log(1./[P(:) P(:)]'),repmat(ylim,length(P),1)','k'); hold off

%
%
%% Move the regression line around by changing slope and intercept
% The red line IS the best-fit across all of the spectrum. But the spectrum is noisy! Move it up and down until you see that the significant peak (loosely defined) really rise above the background.

% You will want to experiment with what are reasonable values, to your eye
offs=0.0;
offi=4.0;
% Construct the new line
yy=polyval([slope+offs intrc+offi],log(fax));
% Plot it, in green this time
hold on; py=plot(log(fax),yy,'g','linew',2); hold off

%
%
%% Zoom in on the interval of interest and add explicit labels
% Find what works to make your point!

% Adjust the zoom
xlim([-4 -1])
ylim([5 25])
% Label explicitly
for index=1:length(P)
    text(log(1/[P(index)])+0.05,22,sprintf('%i h',P(index)))
end
% Take off the red line, it has done its job!
delete(findobj('Color','r'))
% Big up the labels, you'll see they're quite small
xx.FontSize=16;
yy.FontSize=16;
tt.FontSize=16;

%
%% Your job: fjsimonsl04
% Present a Fourier-analysis  of the Hargraves Hall temperature/weather record, broken down any which you want. Such a study could involve the other variables in the weather record, any sort of time-segmenting that you may want to do (do the periodicities vary much week by week, month to month?), and so on (will you want to move the regression line, change the zoom, label more periods?).
%
%% Leave the LiveScript environment!
% Un-anchor the plot by clicking on the encircled arrow in the top right hand of the plot.
% Use the buttons in your figure window to save the plot as a Portable Document Format (*.pdf).
% Alternatively, go to the Command Window, specify a target location (a figure directory), and use fullfile to make a complete path, i.e: type:

% If you see nothing, make sure you click on the figure to make it "active"
% Where will you save this figure? E.g., on your Desktop or somewhere that 
% you can easily find it! I stick it in my Dropbox.
% You may want to do this in the Command Window!
figdir='/Users/fjsimons/Dropbox/FRS124/Labs/L04';
% Issue the print command, which does the same as "save" it
print('-dpdf',fullfile(figdir,'fjsimonsl04-fig.pdf'))

%
%% Leave the MATLAB environment altogether! Time to learn LaTeX!
% We will now teach you to incorporate this figure into a LaTeX document. The key commands are to generate a float (a "floating" figure environment), within which you issue a LaTeX include-graphics command, and the options that you can give it will determine its size, trim the whitespace, generate a label and make a caption. Within LaTeX, this will look something like this here:

% No more MATLAB! Inside your LaTeX document, something of the kind
%\begin{figure}\centering
%\includegraphics[width=0.65\textwidth,trim=2cm 7cm 2cm 6.5cm, clip]
%{/Users/fjsimons/Dropbox/FRS124/Labs/L04/fjsimonsl04-fig}
%\end{figure}