How can I decompose a signal into square waves? I am dealing with signals that are a superposition of different square waves with different amplitudes and phases. Normally, one would decompose a signal into sine waves with help of the Fourier transform, but in this particular case a decomposition into square waves would be much more effective. A Fourier transform would produce a very complicated spectrum, while a square wave decomposition should give just a few clear lines.

I know that such a decomposition is possible. In fact, I could use any periodic function as a basis for the decomposition and this is mentioned in many texts on the subject. But I could never find a formula or an explicit example for a decomposition into a non-sinusoidal basis.

My approach to decompose a signal consisting of the $N$ samples $x_k$, was to use a DFT-like formula
$$ u_k = \sum_{n=0}^{N-1} x_n \, \mathcal{R}_k(n)$$
where $\mathcal{R}_k$ is a real-valued square wave with a frequency $k$ times the base frequency. But this is certainly not complete, since I do not obtain any phase information for the constituent square waves, and I couldn't invert the procedure.

How can I decompose my signals into square waves with well-defined amplitude and phase?
What is described in the question is very near the Discrete Wavelet Transform (DWT) with the use of the Haar Wavelet.

The DWT decomposes a signal into a sum of dilated and translated orthogonal functions that do not necessarily have to be trigonometric. The DWT does not transform a signal from the time domain to a frequency domain but to a scale space where the "time" dimension is preserved. The Haar wavelet is effectively just one period of a square wave and due to its dilation and replication as the transform progresses it would appear as occurring at different frequencies. For more information on the link between decomposition level and frequency please see this link

Another transform that might be of help here, is the Walsh-Hadamard transform which does exactly that, decompose a signal into a sum of square waveforms which are orthogonal (Please note the sequency there as well).

For a brief example that seems to be near what you are after, please see this link

Hope this helps.