Consider the first-order ODE $$ x'(t) = f(t) \cos(\omega t) - \gamma x(t) $$ whose exact solution is $$ x(t) = x(0) e^{-\gamma t} + \int_0^t e^{-\gamma (t-t')} \cos(\omega t') f(t')~dt' $$ For simple forms of $f(t)$, the time-history integral may be performed and one often finds that there is a non-oscillatory term and two oscillating terms, one in-phase with the drive and one out -of-phase. For example if $f(t) = f_0$ is constant, then $$ x(t) = \left[x(0) - \frac{\gamma}{\gamma^2 + \omega^2}\right] e^{-\gamma t} + \frac{\gamma f_0 \cos(\omega t)}{\gamma^2 + \omega^2} + \frac{\omega f_0\sin(\omega t)}{\gamma^2 + \omega^2}\tag{$\star$} $$ The first term does not oscillate, the second term oscillates in phase with the drive, and the third oscillates out of phase. For more complicated $f(t)$, it seems reasonable to expect similar qualitative behavior and to try to extract each type of term with the ansatz $$ x(t)=x_0(t) + \Re[x_1(t)\exp(-i\omega t)] $$ where $x_0$ and $x_1$ are supposed to be non-oscillatory. Putting this in to the original ODE, I obtain $$ x_0'(t) + \gamma x_0(t) + \Re\left[\left(x_1' + \gamma x_1 - i\omega x_1 - f\right) \exp(-i\omega t)\right] = 0 $$ My instinct is to attempt to treat each "Fourier component" independently $$\begin{align} & x'_0(t) + \gamma x_0(t) = 0 \\ & x'_1(t) + \gamma x_1(t) - i\omega x_1(t) = f(t) \end{align}$$ However, this appears not to work as hoped, because the solution for $x_1$ will clearly have an oscillation at the drive frequency $\omega$, due to the term $-i\omega x_1$. For the constant drive case, I have explicitly $$\begin{align} & x_0(t) = x_0(0) e^{-\gamma t} \\ & x_1(t) = \left[x_1(0) - \frac{(\gamma + i\omega)f_0}{\gamma^2+\omega^2} \right] e^{-\gamma t} e^{i\omega t} + \frac{(\gamma + i\omega)f_0}{\gamma^2 + \omega^2} \end{align}$$ In this case, I can actually eliminate the problematic oscillatory term in $x_1$ by choosing to split the initial conditions according to $$ x_0(0) = x(0) - \frac{\gamma f_0}{\gamma^2 + \omega^2}, \quad x_1(0) = \frac{\gamma + i\omega}{\gamma^2 + \omega^2} f_0 $$ In which case I obtain the full solution ($\star$) for $x(t)$, now with the non-oscillatory and oscillatory parts cleanly separated.
So the question I'm interested in is: does this procedure generalize for any (reasonable) $f(t)$? That is, suppose I cannot analytically solve the $x_1$ equation, is there a way to pick out the oscillating and non-oscillating parts? An approximate or asymptotic method would be welcome as well.