I rencently stumbled on a more or less similar question. Actually the Fourier transform of the velocity autocorrelation does not gives you the phonon DOS but rather the phonon population of your system. In other words it gives you normal modes.
The actual calculation of the DOS can be found, for example, in the thesis of Dr <PERSON>, gracefully linked by Prof. <PERSON> in this discussion: https://www.researchgate.net/post/How_do_I_calculate_Phonon_Density_of_states_from_VACF
As Prof. <PERSON> says, the actual calculation is rare in the litterature and it could clarify some question to develop it here.
Unfortunately it is in spanish but I can give here a quick explaination of his calculations.
A. Computation of the phonon DOS
Given the Fourier transform of the velocities from the space of time ($t$) to the space of frequencies ($\omega$):
$$\mathbf{v}n(\omega)=\int{-\infty}^{\infty}\mathbf{v}_n(t)e^{i\omega{}t}dt$$
where the subscript $n$ corresponds to the $n^{th}$ atom and $i$ is the imaginary unit. 346
From here one can get the spectrum of potential kinetic energy:
$$|\mathbf{v}n(\omega)|^2 = \int{-\infty}^{\infty}\int_{-\infty}^{\infty} \mathbf{v}_n(t') \mathbf{v}_n(t) e^{i\omega{}(t-t')}dtdt'$$
In a stationnary oscillatory state, the coordinate $r_n(t)$ can be rewritten as a function of the normal modes of oscillation:
$$ r_{nj}(t) = \sum_s Q_{snj} e^{-i\omega_{sj} t} $$
with the $\omega_s$ the normal frequencies of oscillation and j the direction of the 3D space, $Q_{snj}$ being the mean coordinate of the particle. Hence one can express velocities as the derivative regarding time:
$$ v_{nj}(t) = \sum_s Q_{snj}(-i\omega_{sj}) e^{-i\omega_{sj} t} $$
Using this expression in the second integral, one can deduce:
$$ \sum_{i=1}^N |\mathbf{v}n(\omega_s)|^2 = \sum_s\sum{(n*j)=1}^{3N} \int_{-\infty}^{\infty} |Q_{snj}|^2\omega_s^2 e^{i(\omega+\omega_s)t''}dt'' $$
where $t''=t-t'$.
At thermal equilibrium, given the equipartition of the energy, one has directly $\sum_j |Q_{snj}|^2\omega_{sj}^2=3k_BT$. As the normale modes are supposed to be dominant, the integral of the complex exponential can be reduced to a Dirac distribution such as $\rho(\omega) = \sum_s \delta(\omega+\omega_s)$ is the phonon DOS.
Hence: $$ \sum_{i=1}^N |\mathbf{v}n(\omega_s)|^2 = 3Nk_BT\sum_s\delta(\omega+\omega_s) $$ $$ \downarrow $$ $$ \rho(\omega) = \frac{\sum{i=1}^N |\mathbf{v}_n(\omega_s)|^2}{3Nk_BT} $$
Taking $t'=0$ for the origin of the time, this final expression reduces to $$ \rho(\omega) = \frac{\int_{-\infty}^{\infty}\sum_{i=1}^N <\mathbf{v}_i(t)\cdot\mathbf{v}_i(0)>}{3Nk_BT} $$ linking it to the autocorrelation. If you just Fourier transform the autocorrelation of the velocities, you get the normal modes. You have to had them (i.e integrate) to get the actual DOS. More informations on the assumptions are in the thesis. 37