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?