Time-frequency analysis of non-sinusoidal periodic signals
Let us start with the unsupervised methods...

A first approach would be to compute a spectrogram and factorize it with NMF (Non-negative Matrix Factorization). If you are unfamiliar with this technique, it decomposes a spectrogram into a sum of $k$ constant-spectrum sources, each of these having a time-varying amplitude envelope applied to them. This model perfectly suits your problem and there is a good chance that the columns of the decomposition will be spectra for each  waveform/frequency pair (from which you can estimate a fundamental frequency through spectral sum or by turning it back into an autocorrelation function), and the rows will be the activation signals. There are many implementation out there and it's cheap to just throw your problem at it, there's a good chance it'll just work.

Note that the underlying model of NMF is less restrictive than your signal model...

First, because your activation envelopes are either 0 or 1. I haven't dealt with such constraints in the past, but there's probably a way of defining a penalty measure on the terms of the activation matrix $W$ different from 0 and 1, adding that to the optimized criterion and derive a new set of multiplicative update equations. One constraint not mentioned in your comment is that you probably also assume that the sources are not "blinking" in and out rapidly, and stay active/inactive, for a rather long number of consecutive frames. Constraints can be added to NMF to penalize discontinuous activation envelopes. See Virtanen's paper on continuity constraints for that.

The second specificity of your problem is that you might want to make the assumption that all sources have the same waveform, that is to say, that all sources have the same log-frequency spectrum, modulo a translation. To address this specifity, a recommended technique is to compute a constant-Q spectrogram on your input data (which has the right shift-invariance), and perform a shifted-NMF, with one target source. This will recover the spectrum common to all your sources; and the spectral shifts and temporal activation matrix will provide the equivalent of a spectrogram tailored to your source spectrum. See the synthetic signal in the results section of this paper - this is exactly the signal model you describe. This might require some work, but you could probably tweak this method as well to your more rigid signal model in which the activations are either 0 or 1.

Now on to the supervised methods.

If the set of candidate waveforms is known in advance and small and/or the number of frequencies at which they can appear is small (say you want to transcribe music played on a 80s toy keyboard), you could afford matching pursuit. This will correlate your signal with impulse responses consisting of a windowed signal for all searched waveforms/frequencies pair. Computational costs might rapidly make it a bad choice.

Another supervised method is to use the same kind of multiplicative gradient updates as for NMF, but keeping the matrix $W$ containing the spectral templates "locked". That is to say, you build a matrix W with all candidate spectra you want to probe, and perform your multiplicative update to decompose your spectrogram into a product of an activation matrix H and W. Again, update rules tweakable to include the binary activation constraints. See Bertin's work on piano transcription which uses this technique to factorize spectrograms into a sum of synthetic piano note spectrograms (a good start here). This has also been used for drum transcription by Paulus - using drum samples as a basis.