Infrared and Raman spectroscopies are ubiquitous techniques employed in many experimental laboratories, thanks to their fast and non-destructive nature able to capture materials' features as spectroscopic fingerprints. Nevertheless, these measurements frequently need theoretical and computational support in order to unambiguously decipher and assign complex spectra. Linear-response theory provides an effective way to obtain the higher-order derivatives needed, but its applicability to modern exchange-correlation functionals and pseudopotential formalism remains limited. Here, we devise an automated, open-source, user-friendly approach based on density-functional theory and the electric-enthalpy functional to allow seamless calculation from first principles of infrared absorption and reflectivity, together with zone-center phonons, static dielectric tensor, and Raman spectra. By employing a finite-displacement and finite-field approach, we allow for the use of any functional, as well as an efficient treatment of large low-symmetry structures. Additionally, we propose a simple scheme for efficiently sampling the Brillouin zone at different electric fields. To demonstrate the capabilities of the present approach, we study ferroelectric LiNbO3 crystal as a paradigmatic example, and predict infrared and Raman spectra using various (semi)local, Hubbard corrected, and hybrid functionals. Our results also show how PBE0 and extended Hubbard functionals (PBEsol+U+V) yield for this case the best match in term of peak positions and intensities, respectively.