How can I "smooth" square wave/pulse data to remove/reduce harmonics? First, I have reviewed close-match questions, but can't find something that seems to match what I'm looking for.  Also, it's been 30 years since college, so forgive me for forgetting stuff.. :)

I am analyzing IP Packet data.  I am creating bins (of varying window sizes.. The window seems not to be sufficient to remove/smooth the harmonics) to count the frequency of certain attributes/behaviors within the data.  From this, I am trying to determine fundamental frequencies that activities are occurring.  I believe that the artifacts that I am seeing are harmonics since the data tends to be pulse/bursty, resulting in what amounts to a square wave.  (see below) . In this sample data, as you can see, there is a strong signal occurring every hour, over the course of 47 hours, yet it is almost impossible to "see" because of the harmonics. (at least, that's my guess/intuition)



Is there something that I can do to smooth things out so that the harmonics are suppressed (or largely suppressed) and I only see the fundamental frequencies?  I feel as though I've forgotten something basic.
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);