Higher order nonlinear material mixtures provide a good model to explain the effects of physical–chemical phenomena on hyperspectral remote sensing measurements. Therefore, inverting nonlinear effects starting from the measured spectral values is a very challenging yet fundamental task to provide a thorough and reliable characterization of the materials in a scene. In this paper, this task is achieved by inverting a new model for nonlinear hyperspectral mixtures. Specifically, we show that it is possible to effectively unmix hyperspectral data by assuming a harmonic description of the higher order nonlinear combination of the endmembers. The rationale for this model is that the harmonic analysis is able to understand and quantify effects that cannot be effectively described by classic polynomial combinations. Although the model is nonlinear, unmixing is performed by solving a linear system thanks to the recently proposed polytope decomposition (POD). Experimental results show that inverting this model leads to improved performances with respect to the state of the art in terms of endmember abundance estimation both over synthetic and real datasets.