If you have the Signal Processing Toolbox, I'd use a Savitzky-Golay filter, sgolayfilt, with a third order polynomial and a window size of about 3 months. This will give you something like a sine wave. Then you can subtract that smoothed signal from your original signal to get the hourly deviation from the "average" for that time.
Try this:
fontSize = 16;
s = load('Temperatures.mat')
temps = s.Temp;
subplot(3, 1, 1);
plot(temps, 'b-')
xlabel('Hour', 'FontSize', fontSize)
ylabel('Temperature', 'FontSize', fontSize)
grid on;
title('Actual Temperature', 'FontSize', fontSize)
% Smooth the signal.
windowLength = round(numel(temps) / 60);
% Make sure it's odd
if rem(windowLength, 2) == 0
windowLength = windowLength + 1
end
smoothedTemps = sgolayfilt(temps, 2, windowLength);
subplot(3, 1, 2);
plot(smoothedTemps, 'r-');
xlabel('Hour', 'FontSize', fontSize)
ylabel('Temperature', 'FontSize', fontSize)
title('Seasonal Average Temperature', 'FontSize', fontSize)
grid on;
% Get the deviation from average.
deviationtemps = temps - smoothedTemps;
subplot(3, 1, 3);
plot(deviationtemps, 'b-');
xlabel('Hour', 'FontSize', fontSize)
ylabel('Temperature', 'FontSize', fontSize)
grid on;
title('Detrended Temperature', 'FontSize', fontSize)
If you want, you can get the average of all the years and synthesize a smoothed signal from that. Or you can use interp1 to fill in gaps in your original data.
