How can I "smooth" square wave/pulse data to remove/reduce harmonics?
If you are interested in finding the periods, have a look at autocorrelation lag.
Maybe try some xcorr options too.
MATLAB example code:

nh = 47;  % 47 hours.
% Pulse period of 1 hr and make sequence length not a multiple of period
n = 3600*nh + 273; % 1 sec sampling rate
x = zeros(1,n);
% Create a pulse train with varying heights and some noise
x(777:3600:end) = 8 + 2*rand(1,nh);
x = x + rand(1,n);
% Look at lags up to 4 hours
xcorr(x, 4*3600, 'unbiased'); 


Zero lag is in the middle of the figure.


Note that the consistent heights show reasonably consistent correlation at integer multiples of the base period.
If you want to tolerate timing jitter in the period, you can filter the data with a lowpass FIR filter first, which will smear the peaks in the data and, of course, smear the peaks in the autocorrelation., e.g., 

b = hamming(15);  % Smooth, simple, positive, symmetric, integer delay
x = filter(b,1,x);