Context. X-ray pulsars are binary systems which consist of a neutron star in orbit with a mass donor (companion). In these systems the neutron accretes matter from the companion star, which creates accretion columns or hot spots on the neutron star surface and gives rise to pulsations in the X-ray light curve. The pulse profiles carry information about the accretion and magnetic field geometry.Here we present a study and classification of energy-resolved pulse profiles of a sample of X-ray pulsars, focusing on high-mass X-ray binaries. Aims. Our goal is to perform a classification of X-ray pulsars based on their observed pulse profiles and look for correlations between this classification and their principle physical observables. The analysis pipeline is available online. Methods. We analysed the pulse profiles of a sample of X-ray pulsars using data obtained with the X-ray Multi-Mirror Mission (XMM-Newton) and the Nuclear Spectroscopic Telescope Array (NuSTAR ). We fitted the energy-resolved pulse profiles with a Fourier series of up to five harmonics. We then used the energy evolution of the different Fourier components to classify the pulse profiles into groups. We investigated relationships between the pulse profile properties and other observables of the systems (e.g. orbital period, magnetic field strength, and luminosity) to study the extreme physics of these systems.Results. The sources were divided into three groups using a classification based on the shape, the dominance of the fitted Fourier harmonics, and their respective evolution with energy. We do not find a conclusive correlation between the pulse profile shapes or groups and other parameters of the systems. However, a weak trend was found when comparing our classification to the sources' locations in the spin period-orbital period diagram. Further studies are required to confirm this trend. Conclusions. Despite the large variety of pulse profiles of the X-ray pulsars, we found that with our approach clear categories emerge which we use to classify their behaviour as a function of energy. As we do not find a clear relationship between our classification scheme and other parameters, like the luminosity, the magnetic field strength, or the orbital and spin periods, we conclude that X-ray pulse profiles are influenced by other hidden variables.