A new method in the frequency domain for the identification of nonlinear vibrating structures is described, by adopting the perspective of nonlinearities as internal feedback forces. The technique is based on a polynomial expansion representation of the frequency response function of the underlying linear system, relying on a z-domain formulation. A least squares approach is adopted to take into account the information of all the frequency response functions but, when large data sets are used, the solution of the resulting system of algebraic linear equations can be a difficult task. A procedure to drastically reduce the matrix dimensions and consequently the computational cost – which largely depends on the number of spectral lines – is adopted, leading to a compact and well conditioned problem. The robustness and numerical performances of the method are demonstrated by its implementation on simulated data from single and two degree of freedom systems with typical nonlinear characteristics.