Context
Consider an $N$-level Hamiltonian with energies $\omega_1...\omega_N$ with coupling drives at frequencies $f_{i,j}$ which couple the $i$ and $j$-th levels (not necessarily resonantly, so $f_{i,j} \neq \omega_j - \omega_i$).
This Hamiltonian looks something like
$$H = \sum_{i}^{N} \omega_i \vert i\rangle \langle i \vert + \sum_{i,j} \Omega_{i,j}\cos(f_{i,j} t) \vert i\rangle \langle j \vert + h.c.$$ Or as a matrix
$$H= \begin{bmatrix} \omega_1 & \Omega_{1,2}\cos(f_{1,2} t) & \Omega_{1,3}\cos(f_{1,3} t) & \dots\\ \Omega^*_{1,2}\cos(f_{1,2} t)& \omega_2& \Omega_{2,3}\cos(f_{2,3} t) & \dots \\ \Omega_{1,3}^*\cos(f_{1,2} t) & \Omega^*_{2,3}\cos(f_{2,3} t)&\omega_3 & \dots\\ \vdots & \vdots &\vdots & \ddots \end{bmatrix}$$
The rotating wave approximate to $H$, denoted $H_{RWA}$ would involve first replacing the $\cos(f_{i,j} t)$ terms with $e^{i f_{i,j} t}$
$$H\approx \tilde{H}=\begin{bmatrix} \omega_1 & \Omega_{1,2}e^{if_{1,2} t} & \Omega_{1,3}e^{if_{1,3} t}& \dots\\ \Omega^*_{1,2}e^{-if_{1,2} t}& \omega_2& \Omega_{2,3}e^{if_{2,3} t} & \dots \\ \Omega_{1,3}^*e^{-if_{1,2} t} & \Omega^*_{2,3}e^{-if_{2,3} t}&\omega_3 & \dots\\ \vdots & \vdots &\vdots & \ddots \end{bmatrix}$$
The second step is a matter of convenience which is to come up with an appropriate unitary matrix $U$ to get rid of the explicit time dependent with
$$H_{RWA}= U^{\dagger} \tilde{H} U + iU^{\dagger}\partial_t U$$
Question
Is there an algorithm for going through the frequencies appearing in the Hamiltonian $H$ and coming up with the appropriate rotating frame transform $U$ and obtain $H_{RWA}$ from $H$?
It is easy enough to get $U$ with the RWA when you have two levels, but it is unclear to me how to do it systematically for larger $N$.
