It is a well-known result that when two signals with the same period $T$ get multiplied in the time domain, the resulting signal's Fourier coefficients are given by the discrete convolution of individual Fourier coefficients. $$c_k = \sum_{n=-\infty}^{+\infty} a_nb_{k-n}$$
where $c_k$ is the resulting signal's Fourier coefficients and $a_n$ and $b_n$ are the constituent signal, Fourier coefficients. What I am interested in knowing is if the same is true for two signals with different frequencies.
To start off, the two frequencies should at least be rational multiples as explained here. So, if we assume $\omega_x = p\omega_0$ and $\omega_y = q\omega_0$ and follow the steps for inspecting the nature of the resulting signal's fourier coefficients by,
$$c_k = \frac{1}{T}\int_0^T \sum_{n=-\infty}^{+\infty} a_{n} e^{jpn\omega_0t}\sum_{l=-\infty}^{+\infty}b_{l}e^{jlq\omega_0t} e^{-jk\omega_0t}dt$$ Then by following the procedure shown here, we arrive at,
$$c_k = \frac{1}{T} \int_0^T \sum_{n=-\infty}^{+\infty}a_{n}\sum_{l=-\infty}^{+\infty} b_{l}\delta(k-(np+lq))$$
Normally, we would apply sifting property here if $p =1$ and $q=1$ and proceed to get, $$c_{k}= \frac{1}{T} \sum _{n}a_{n} b_{k-n} = \mathbf{a}*\mathbf{b}$$
However, I am unsure how to proceed with the case where we have $\delta(k-(np+lq))$. How do we go about applying time sifting property to the above delta function to arrive at a definitive result (if at all possible)?