We develop an ab initio approach for large-scale Raman scattering spectra simulations based on time-dependent density functional theory in conjunction with the plane-wave pseudopotential method at the Gamma point. A Lagrangian functional is introduced to analytically compute the first-order derivatives of the frequency dependent polarizability with respect to nuclear coordinates. The computational effort of Raman intensities required by our method is reduced by one power of system size compared to that required for phonon frequency calculations. The method is validated for several molecular and solid systems including CH4, C2H2, C6H6, C60, bulk Si, and boron doped Si crystal and exhibits excellent agreements with the experimental Raman spectra. We show that the method yields a computational scaling of Ne2, with Ne ranging from 32 to 4000 electrons, opening doors for many large-scale Raman spectra computations that are beyond the reach of previous approaches.