Decompose a time series into superposition of step functions? Background
I have time series data comprising hourly observations of a sensor's readings over a period of almost a year. The sensor records an environment whose baseline measurements should have daily cyclicality, and which are generated by a sum of several 'steppy' (i.e. on, off) processes. However, there is a lot of other stuff going on in the environment, and the sensor captures that as well. The 'other stuff' often has a much larger amplitude than the underlying cyclical data.
A final caveat - the 'other stuff' may actually have roughly daily cyclicality, but if so its period is extremely noisy - i.e. the 'other stuff' (say) starts at 1200 on average, and ends at 1500 on average, but the start and end points day to day vary widely.
Question


*

*How can I uncover the actual cyclical step processes underlying this data? (main question)

*Can I write a sum-of-squared-error based optimisation that will try and represent this data as a superposition of step functions? I'm thinking something like multiple linear regression except the IDVs encode step functions with different amplitudes and starting times? (secondary question) 


Some approaches that either havent worked, or that I dont think will work


*

*Fourier analysis. I've run DFT, and tried saving either the largest amplitude frequencies, or the ones whose periods are between 0.8 and 1 days. However, this doesn't successfully filter out the 'other stuff' - in fact, because the base cyclical process is steppy, retaining only a few frequencies cant really capture this adequately, and because the SNR ratio is low, the reconstructed series can be quite close to the 'single-day-average' of the noise process

*See 2. above in Questions. I don't think this will work. Two
appoaches: 1) pretend I'm doing simple linear regression, with one
IDVs that encode every possible step function (e.g if a day was
represented by only 4 time points, {1, 1, 0, 0}, {0, 0, 1, 1}, ... etc
(likely to be overparameterized) 2) or write a sum of squares cost function where my estimator is a sum over S step functions, parameterized by an amplitude and a step time. I'm not sure that if I optimise my cost function against such an esimator that it will converge to a global minimum, especially if the number of steps, S, is allowed to vary.
Would really appreciate some insight here!
Thanks
 A: If you have two signals with the same frequency, it's going to be difficult to separate them out. One idea is separation of the phases, maybe your environment "noise" is phase-shifted to the signal. 
Note, that I'm using quotes for noise, because it really is not. The noise must be white, while yours clearly is not. So, technically, it's a signal, which you want to separate out.
Try fourier transform, but instead of, as usually, looking at spectral densities, pay attention to phases.
The second idea is to use some sort of wavelet transforms. If you use straight FFT(DFT) there will be many overtones from both the signal and the "other stuff", which may make it difficult to separate out. Maybe wavelets would fit better in this case.
UPDATE: In any case if your signal and other stuff have the same base frequency, it's a huge issue. Think of an extreme case: two sines with the same frequency and different phases: $\sin(wt)+\sin(wt+\phi)=2\cos(\phi/2)\sin(wt+\phi/2)$. This is just precisely one sine wave, FFT will not be able to detect two waves.
A: Fourier transforms and other periodograms are mostly good if you do not know the underlying period of your data. (Moreover 24 data points per cycle with noise will make it very difficult to gather the subharmonics etc. needed to reconstruct the underlying periodic signal.)
The main advantage that you have is that you know the periodicity of the data or have some hypothesis about it that you want to test. Thus I would suggest two much simpler approaches:


*

*Take averages of all the data for a given time of the day and look at the resulting time series, which represents the average daily cycle. If your assumptions are correct, the effect of the noise should be reduced by the averaging and your actual steps should remain.

*Take a look at all pairs of adjacent times of the day and count the number of increases (or decreases) that exceed a given threshold. If your assumptions are correct and your noise level is not too high, there should be a threshold for which your count should be significantly higher for some times of the day.
If both of the above approaches fail, your hypothesis (as I understood it) is likely not correct.
A: Have you tried using periodogram methods, such as Welch's method? If the 'other stuff' is variable in frequency over time, where as the desired signal stays constant in frequency over time, then there are better estimates of the spectrum than DFT on the entire dataset. Welch's method takes several DFTs on a sliding window of your data, and averages them together.
I don't know enough about your data to answer your second question. Perhaps you should post an example plot? If you have several step functions then the nonlinear least squares model will take forever to solve unless there is another relationship in your data that I am not aware of.
