A technique is reported for generating the quantitative area under selected peaks within a 31P NMR spectrum. A numerical iterative method generates the fitted curve so as to minimize the RMS deviation between the fit and the experimental data. The curve is constructed from elemental grains of spectral density ('spexels'), each of which represents an elemental Lorentzian distribution, where the centre frequency and line width of the spexel may be varied within predetermined limits. This provides a fit that in principle is not restricted to a Lorentzian model. The method allows peak areas to be estimated, including the case of overlapping peaks. The method has been tested using simulated spectra containing six overlapping spectral lines each of known amplitude (ranging from 367 to 661 mV) and area; together with additional Gaussian noise with a standard deviation ranging from 30 to 646 mV. The results of fitting both unfiltered and filtered spectra were compared. The variation of quality of fit with spectral noise and filtering has been evaluated. In all cases, the method fitted the peak amplitudes to within 1% of the simulated value. The optimization processes provide an excellent non-linear spectrum filtering algorithm. Provided the noise in the spectra did not exceed ca 400 mV, prefiltered data could be adequately fitted.