We introduce a fast Fourier spectral method to compute the linearized Boltzmann collision operators of the Boltzmann equation for all variable hard-sphere gases. The linearization is performed about any local Maxwellian, so that the resulting model can preserve the Euler and Navier-Stokes limits. While the state-of-the-art method provides a computational complexity O(M N^4 log N), with N being the number of modes in each direction and M being the number of quadrature points on a hemisphere, our method removes the factor M, reducing the computational cost to O(N^4 log N). The method is applied to a numerical solver for the steady-state Boltzmann equation with the quadratic collision operator, using a modified Newton's method where the Jacobian of the collision term is approximated by a linearization of the collision operator about the local Maxwellian. Our numerical experiments show that the new method maintains the fast convergence of Newton's method, while the computational cost for solving the linear system can be significantly reduced.