Fast multidimensional NMR is important in chemical shift assignment and for studying structures of large proteins. We present the first method which takes advantage of the sparsity of the wavelet representation of the NMR spectra and reconstructs the spectra from partial random measurements of its free induction decay (FID) by solving the following optimization problem: min ‖ x ‖ 1 subject to ‖ y − S F T W T x ‖ 2 ≤ ε , where y is a given n × 1 observation vector, S a random sampling operator, F denotes the Fourier transform, and W an orthogonal 2D wavelet transform. The matrix A = S F T W T is a given n × p matrix such that n < p . This problem can be solved by general-purpose solvers; however, these can be prohibitively expensive in large-scale applications. In the settings of interest, the underlying solution is sparse with a few nonzeros. We show here that for large practical systems, a good approximation to the sparsest solution is obtained by iterative thresholding algorithms running much more rapidly than general solvers. We demonstrate the applicability of our approach to fast multidimensional NMR spectroscopy. Our main practical result estimates a four-fold reduction in sampling and experiment time without loss of resolution while maintaining sensitivity for a wide range of existing settings. Our results maintain the quality of the peak list of the reconstructed signal which is the key deliverable used in protein structure determination.