The article is devoted to methods of discretization of energy characteristics of two-dimensional random signals when simulating random signals using the original harmonic method, which is a generalization of the well-known algorithm proposed by V. S. Pugachev for the two-dimensional case. Requirements imposed on the sampling method are aimed at reducing the computational complexity of the simulation method and increasing its flexibility thanks to removing restrictions on the form of autocorrelation functions and spectral energy density functions. The use of the simulation error as a criterion for quality assessment is proposed. The discretization method is considered for signals given both on unlimited definition intervals and on limited ones. The article demonstrates results of the software system implementation in which the original simulation method is realized using the described sampling methods in both cases. The proposed technique is shown to be robust and efficient, with the results obtained being of independent scientific and technical value and showing promise for developing new effective spectral techniques of simulating signals for the use in intelligent decision support systems.