We have developed a formulation of density-functional perturbation theory for the calculation of vibrational frequencies in molecules and solids, which uses numerical atomic orbitals as a basis set for the electronic states. The (harmonic) dynamical matrix is extracted directly from the first-order change in the density matrix with respect to infinitesimal atomic displacements from the equilibrium configuration. We have applied this method to study the vibrational properties of a number of hydrogen-related complexes and light impurities in silicon. The diagonalization of the dynamical matrix provides the vibrational modes and frequencies, including the local vibrational modes (LVM’s) associated with the defects. In addition to tests on simple molecules, results for interstitial hydrogen, hydrogen dimers, vacancy-hydrogen and self-interstitial-hydrogen complexes, the boron-hydrogen pair, substitutional C, and several O-related defects in c-Si, are presented. The average error relative to experiment for the ∼60 predicted LVM’s is about 2% with most highly harmonic modes being extremely close and the more anharmonic ones within 5–6 % of the measured values.