One can characterize long-memory Gaussian stationary processes by sparse covariance of wavelet coefficients. However, these covariance moments do not capture statistical dependencies of the Fourier coefficients of a non-Gaussian process. Such dependencies manifest in the phase of complex-valued wavelet coefficients. Wavelet phase harmonics multiply the phase of wavelet coefficients with integers. By computing their covariance, we are able to capture key non-Gaussianities of a stationary process including high-order structure functions in Turbulence analysis. This is validated by sampling a maximum-entropy model conditioned on a sparse subset of this covariance. Given a finite number of data samples of a non-Gaussian stationary process such as texture and turbulent fluids, we demonstrate the accuracy of the maximum-entropy models through a fundamental approximation and estimation trade-off using the method of moments.