Interpreting Raman and IR vibrational spectra in complex organic molecules lacking symmetries poses a formidable challenge. In this study, we propose an innovative approach for simulating vibrational spectra and attributing observed peaks to molecular motions, even when highly anharmonic, without the need for computationally expensive ab initio calculations. Our approach stems from the time-dependent stochastic self-consistent harmonic approximation to capture quantum nuclear fluctuations in atom dynamics while describing interatomic interaction through state-of-the-art reactive machine-learning force fields. Finally, we employ an isotropic charge model and a bond capacitor model trained on ab initio data to predict the intensity of IR and Raman signals.