A global stability analysis of a Falkner-Skan-Cooke boundary layer with distributed three-dimensional surface roughness is per- formed using high-order direct numerical simulations. Computations have been performed for different sizes of the roughness elements, and a time-stepping method has been used to find the instability modes. The study shows that a critical roughness height beyond which a global instability is excited does exist. Furthermore, the origins of this instability is examined by means of an energy analysis, which reveals the production and dissipation terms responsible for the instability, as well as the region in space where the instability originates.