We propose a step-by-step experimental procedure for characterization of the nonlinear contact stiffness on surfaces using contact-mode atomic force microscopy. Our approach directly estimates the first-, second-, and third-order coefficients of the contact stiffness. It neither uses nor requires the underlying assumptions of the Hertzian contact theory. We use a primary resonance excitation of the probe to estimate the linear coefficient of the contact stiffness. We use the method of multiple scales to obtain closed-form expressions approximating the response of the probe to a subharmonic resonance excitation of order one-half. We utilize these expressions and higher-order spectral measurements to independently estimate the quadratic and cubic coefficients of the contact stiffness.