A technique based on polarization modulation spectroscopy (PMS) has been developed to determine quantitatively the number of fluorophores in nanoparticles at the single-molecule level. The technique involves rotation of the polarization of the excitation laser on a millisecond time scale, leading to fluorescence intensity modulation. By taking account of the heterogeneous orientation among the dipoles of the fluorophores and simulating the modulation depth distribution with Monte Carlo calculations, we show that it is possible to deduce the ensemble average and number distribution of the fluorophores. We apply the technique to fluorescent nanodiamonds (FNDs) containing multiple nitrogen vacancy (NV) centers. Comparing the experimental and simulated modulation depth distributions of 11 nm FNDs, we deduce an average number of <N> = 3, which is in good agreement with independent photon correlation measurements. The method is general, rapid, and applicable to other nanoparticles, polymers, and molecular complexes containing multiple and randomly orientated fluorophores as well.