We present a practical, ab initio time-dependent method using Hagedorn wavepackets to efficiently simulate single vibronic level (SVL) fluorescence spectra of polyatomic molecules from arbitrary initial vibrational levels. We apply the method to compute SVL spectra of anthracene by performing wavepacket dynamics on a 66-dimensional harmonic potential energy surface constructed from density functional theory calculations. The Hagedorn approach captures both mode distortion (frequency changes) and mode mixing (Duschinsky rotation) within the harmonic approximation. We not only reproduce the previously reported simulation results for singly excited 121 and 1̅1̅1 levels but are also able to compute SVL spectra from multiply excited levels in good agreement with experiments. Notably, all spectra were obtained from the same wavepacket trajectory without any additional propagation beyond what is required for the emission spectrum from the ground vibrational level of the electronically excited state.