The Fermi-Dirac distribution gives the probability of an energy level being occupied. You have to combine it with the density of states $g(\epsilon)$ at that energy to get the physical properties you could be interested in, for example, the internal energy $U$ can be computed as \begin{equation} U = \int \epsilon\, f_{FD}(\epsilon)\, g(\epsilon)\,\mathrm{d}\epsilon \end{equation}

The expression of the density of states depends on the system you are considering.

In your formula, the prefactor is indeed due to the density of states of a gas of free electrons. You can find the detailed calculations in many solid-state textbooks (a personal favorite is Ashcroft-Mermin).

The calculation is something along the following lines:

1. Consider a sample of volume $V$ in real space. The number of free-electrons levels in the momentum space, or $\mathbf{k}$-space, about a point $\mathbf{k}$ of volume $\mathrm{d}\mathbf{k}$ for this sample is \begin{equation} \mathrm{d}N = 2\frac{V}{8\pi^3}\,\mathrm{d}\mathbf{k}, \end{equation} where the 2 is due to the twofold spin degeneracy of spin-1/2 particles.
2. The free-electron energy density in the $\mathbf{k}$-space, can be written as \begin{equation} g(\epsilon(\mathbf{k})) =\frac{1}{V}\frac{\mathrm{d}N}{\mathrm{d}\mathbf{k}}= 2 \frac{1}{8 \pi^3} \end{equation}
3. Of course, for the free electron gas, the energy corresponds to the kinetic energy \begin{equation} \epsilon(\mathbf{k}) = \frac{\hbar^2\mathbf{k}^2}{2m} \end{equation}
4. If we multiply the density of states $g(\epsilon(\mathbf{k}))$ with the Fermi-Dirac distribution (= the probability of each level to be occupied), we obtain the density of electrons with momentum in the volume element $\mathrm{d}\mathbf{k}$ per unit volume of real space \begin{equation} F(\mathbf{k}) = g(\epsilon(\mathbf{k})) f_{FD}(\epsilon(\mathbf{k})) = 2 \frac{1}{8 \pi^3} \frac{1}{\exp\left(\frac{\frac{\hbar^2\mathbf{k}^2}{2m} - \mu}{k_BT}\right) + 1} \end{equation}
5. As the number of particles must be conserved, we must have that the number of electrons in an element of volume $\mathrm{d}\mathbf{v}$ about $\mathbf{v}$ is the same as the number of electrons in an element of volume $\mathrm{d}\mathbf{k}$ about $\mathbf{k}$: \begin{equation} F(\mathbf{k})\,\mathrm{d}\mathbf{k} = F(\mathbf{v})\,\mathrm{d}\mathbf{v}, \end{equation} where the velocity is given by $\mathbf{v}=\frac{\hbar \mathbf{k}}{m}$. In cartesian coordinates the above becomes \begin{equation} F(\mathbf{k})\,\mathrm{d}k_x\,\mathrm{d}k_y\,\mathrm{d}k_z = F(\mathbf{v})\,\mathrm{d}v_x\,\mathrm{d}v_y\,\mathrm{d}v_z \end{equation}

From here you can easily find that \begin{equation} F(\mathbf{v}) = F(\mathbf{k})\frac{\mathrm{d}k_x}{\mathrm{d}v_x}\frac{\mathrm{d}k_y}{\mathrm{d}v_y}\frac{\mathrm{d}k_z}{\mathrm{d}v_z} = F(\mathbf{k}) \left(\frac{m}{\hbar}\right)^3 \end{equation}

The Fermi-Dirac distribution gives the probability of an energy level being occupied. You have to combine it with the density of states $g(\epsilon)$ at that energy to get the physical properties you could be interested in, for example, the internal energy $U$ can be computed as \begin{equation} U = \int \epsilon\, f_{FD}(\epsilon)\, g(\epsilon)\,\mathrm{d}\epsilon \end{equation}

The expression of the density of states depends on the system you are considering.

In your formula, the prefactor is indeed due to the density of states of a gas of free electrons. You can find the detailed calculations in many solid-state textbooks (a personal favorite is Ashcroft-Mermin).

The calculation is something along the following lines:

1. Consider a sample of volume $V$ in real space. The number of free-electrons levels in the momentum space, or $\mathbf{k}$-space, about a point $\mathbf{k}$ of volume $\mathrm{d}\mathbf{k}$ for this sample is \begin{equation} \mathrm{d}N = 2\frac{V}{8\pi^3}\,\mathrm{d}\mathbf{k}, \end{equation} where the 2 is due to the twofold spin degeneracy of spin-1/2 particles.
2. The density of allowed free-electron states in the $\mathbf{k}$-space per unit volume, can be written as \begin{equation} g(\epsilon(\mathbf{k})) =\frac{1}{V}\frac{\mathrm{d}N}{\mathrm{d}\mathbf{k}}= 2 \frac{1}{8 \pi^3} \end{equation}
3. Of course, for the free electron gas, the energy corresponds to the kinetic energy \begin{equation} \epsilon(\mathbf{k}) = \frac{\hbar^2\mathbf{k}^2}{2m} \end{equation}
4. If we multiply the density of states $g(\epsilon(\mathbf{k}))$ with the Fermi-Dirac distribution (= the probability of each level to be occupied), we obtain the density of electrons with momentum in the volume element $\mathrm{d}\mathbf{k}$ per unit volume of real space \begin{equation} F(\mathbf{k}) = g(\epsilon(\mathbf{k})) f_{FD}(\epsilon(\mathbf{k})) = 2 \frac{1}{8 \pi^3} \frac{1}{\exp\left(\frac{\frac{\hbar^2\mathbf{k}^2}{2m} - \mu}{k_BT}\right) + 1} \end{equation}
5. As the number of particles must be conserved, we must have that the number of electrons in an element of volume $\mathrm{d}\mathbf{v}$ about $\mathbf{v}$ is the same as the number of electrons in an element of volume $\mathrm{d}\mathbf{k}$ about $\mathbf{k}$: \begin{equation} F(\mathbf{k})\,\mathrm{d}\mathbf{k} = F(\mathbf{v})\,\mathrm{d}\mathbf{v}, \end{equation} where the velocity is given by $\mathbf{v}=\frac{\hbar \mathbf{k}}{m}$. In cartesian coordinates the above becomes \begin{equation} F(\mathbf{k})\,\mathrm{d}k_x\,\mathrm{d}k_y\,\mathrm{d}k_z = F(\mathbf{v})\,\mathrm{d}v_x\,\mathrm{d}v_y\,\mathrm{d}v_z \end{equation}

From here you can easily find that \begin{equation} F(\mathbf{v}) = F(\mathbf{k})\frac{\mathrm{d}k_x}{\mathrm{d}v_x}\frac{\mathrm{d}k_y}{\mathrm{d}v_y}\frac{\mathrm{d}k_z}{\mathrm{d}v_z} = F(\mathbf{k}) \left(\frac{m}{\hbar}\right)^3 \end{equation}