In single vibronic level (SVL) fluorescence experiments, the electronically excited initial state is also excited in one or several vibrational modes. Because computing such spectra by evaluating all contributing Franck–Condon factors becomes impractical (and unnecessary) in large systems, here we propose a time-dependent approach based on Hagedorn wavepacket dynamics. We use Hagedorn functions—products of a Gaussian and carefully generated polynomials—to represent SVL initial states because in systems whose potential is at most quadratic, Hagedorn functions are exact solutions to the time-dependent Schrödinger equation and can be propagated with the same equations of motion as a simple Gaussian wavepacket. Having developed an efficient recursive algorithm to compute the overlaps between two Hagedorn wavepackets, we can now evaluate emission spectra from arbitrary vibronic levels using a single trajectory. We validate the method in two-dimensional global harmonic models by comparing it with quantum split-operator calculations. In addition, we study the effects of displacement, distortion (squeezing), and Duschinsky rotation on SVL fluorescence spectra. Finally, we demonstrate the applicability of the Hagedorn approach to high-dimensional systems on a displaced, distorted, and Duschinsky-rotated harmonic model with 100 degrees of freedom.