Yes, however there are some details to iron out. Essentially, you have the discrete analogue of Fourier transform:
$$
\phi(k)=\sum_{n=-\infty}^{+\infty}\psi_ne^{in k} \\
\psi_n=\int_0^{2\pi}\frac{dk}{2\pi}\phi(k)e^{ink}
$$
You therefore get:
$$
\psi_n(t)=\int_0^{2\pi}\frac{dk}{2\pi}(A_+(k)e^{ink-i\omega(k)t}+A_-(k)e^{ink+i\omega(k)t})
$$
Note that per mode you have two coefficients, which is normal since the equation is second order in time. The coefficients $A_\pm(k)$ can be determined by transforming the initial conditions:
$$
\psi_n(0)=\int_0^{2\pi}\frac{dk}{2\pi}B(k)e^{ink} \\
\dot \psi_n(0)=\int_0^{2\pi}\frac{dk}{2\pi}C(k)e^{ink} \\
$$
Identifying, you get:
$$
A_\pm(k) = \frac{1}{2}\left(B(k)\pm\frac{i}{\omega(k)}C(k)\right)
$$
Btw using the dispersion relation, you can integrate out $k$ to obtain the propagators for initial field and initial derivative.
Hope this helps.
Answer to comment
You can make the dependence of the field in the initial conditions more explicit. You can rewrite it as:
$$
\psi_n(t)=\sum_{m\in\mathbb Z}\psi_m(0)K_{n-m}^0(t)+\dot\psi_m(0)K_{n-m}^1(t)
$$
with $K^0,K^1$ being the propagators or kernels. You don’t need to go to Fourier space if you have their expressions.
Note that computationally, it is faster to use the previous method since Fourier transforming is faster than calculating a convolution. However, the kernels are often useful to understand the large time behavior which is not obvious from the Fourier method.
In our case, using the previous formulas:
$$
K^0_n(t)=\int_0^{2\pi}\frac{dk}{2\pi}e^{ink}\cos(\omega(k)t) \\
K^1_n(t)=\int_0^{2\pi}\frac{dk}{2\pi}e^{ink}\frac{\sin(\omega(k)t)}{\omega(k)}
$$
The large $t$ asymptotics can be recovered using the saddle-point method. However, in your case, I don’t know whether a closed form exists. In the continuum limit, you do get a simple expression though.
