There is a good reason for these definitions, which becomes clearer when you look at the general form for moments of standardised random variables. To answer this question, first consider the general form of the $k$th standardised central moment:$^\dagger$
$$\phi_k = \mathbb{E} \Bigg[ \Bigg( \frac{X - \mathbb{E}[X]}{\mathbb{S}[X]} \Bigg)^k \text{ } \Bigg].$$
The first two standardised central moments are the values $\phi_1=0$ and $\phi_2=1$, which hold for all distributions for which the above quantity is well-defined. Hence, we can consider the non-trivial standardised central moments that occur for values $k \geqslant 3$. To facilitate our analysis we define:
$$\begin{equation} \begin{aligned}
\phi_k^+
&= \mathbb{E} \Bigg[ \Bigg| \frac{X - \mathbb{E}[X]}{\mathbb{S}[X]} \Bigg|^k \text{ } \Bigg| X > \mathbb{E}[X] \Bigg] \cdot \mathbb{P}(X > \mathbb{E}[X]), \\[8pt]
\phi_k^- &= \mathbb{E} \Bigg[ \Bigg| \frac{X - \mathbb{E}[X]}{\mathbb{S}[X]} \Bigg|^k \text{ } \Bigg| X < \mathbb{E}[X] \Bigg] \cdot \mathbb{P}(X < \mathbb{E}[X]).
\end{aligned} \end{equation}$$
These are non-negative quantities that give the $k$th absolute power of the standardised random variable conditional on it being above or below its expected value. We will now decompose the standardised central moment into these parts.
Odd values of $k$ measure the skew in the tails: For any odd value of $k \geqslant 3$ we have an odd power in the moment equation and so we can write the standardised central moment as $\phi_k = \phi_k^+ - \phi_k^-$. From this form we see that the standardised central moment gives us the difference between the $k$th absolute power of the standardised random variable, conditional on it being above or below its mean respectively.
Thus, for any odd power $k \geqslant 3$ we will get a measure that gives positive values if the expected absolute power of the standardised random variable is higher for values above the mean than for values below the mean, and gives negative values if the expected absolute power is lower for values above the mean than for values below the mean. Any of these quantities could reasonably be regarded as a measure of a type of "skewness", with higher powers giving greater relative weight to values that are far from the mean.
Since this phenomenon occurs for every odd power $k \geqslant 3$, the natural choice for an archetypal measure of "skewness" is to define $\phi_3$ as the skewness. (The higher-order odd moments $k=5,7,9,...$ are sometimes called measures of "hyperskewness".)This is a lower standardised central moment than the higher odd powers, and it is natural to explore lower-order moments before consideration of higher-order moments. In statistics we have adopted the convention of referring to this standardised central moment as the skewness, since it is the lowest standardised central moment that measures this aspect of the distribution. (The higher odd powers also measure types of skewness, but with greater and greater emphasis on values far from the mean; these are sometimes called measures of "hyperskewness".)
Even values of $k$ measure fatness of tails: For any even value of $k \geqslant 3$ we have an even power in the moment equation and so we can write the standardised central moment as $\phi_k = \phi_k^+ + \phi_k^-$. From this form we see that the standardised central moment gives us the sum of the $k$th absolute power of the standardised random variable, conditional on it being above or below its mean respectively.
Thus, for any even power $k \geqslant 3$ we will get a measure that gives non-negative values, with higher values occurring if the tails of the distribution of the standardised random variable are fatter. Note that this is a result with respect to the standardised random variable, and so a change in scale (changing the variance) has no effect on this measure. Rather, it is effectively a measure of the fatness of the tails, after standardising for the variance of the distribution. Any of these quantities could reasonably be regarded as a measure of a type of "kurtosis", with higher powers giving greater relative weight to values that are far from the mean.
Since this phenomenon occurs for every even power $k \geqslant 3$, the natural choice for an archetypal measure of kurtosis is to define $\phi_4$ as the kurtosis. This is a lower standardised central moment than the higher even powers, and it is natural to explore lower-order moments before consideration of higher-order moments. In statistics we have adopted the convention of referring to this standardised central moment as the "kurtosis", since it is the lowest standardised central moment that measures this aspect of the distribution. (The higher even powers also measure types of kurtosis, but with greater and greater emphasis on values far from the mean; these are sometimes called measures of "hyperkurtosis".)
$^\dagger$ This equation is well defined for any distribution whose first two moments exist, and which has non-zero variance. We will assume that the distribution of interest falls in this class for the rest of the analysis.
Best Answer
I do not know about any exact matching technique. If an approximated technique could work for you (meaning getting a density which approximately matches m given moments), then you could consider using an orthogonal polynomial series approach. You would choose a polynomial basis (Laguerre, Hermite, etc) depending on the range of your data. I describe below a technique that I have used in Arbel et al. (1) for a compactly supported distribution (see details in Section 3).
In order to set the notation, let us consider a generic continuous random variable $X$ on $[0,1]$, and denote by $f$ its density (to be approximated), and its raw moments by $\gamma_r=\mathbb{E}\big[X^r\big]$, with $r\in\mathbb{N}$. Denote the basis of Jacobi polynomials by $$G_i(s) = \sum_{r=0}^i G_{i,r}s^r,\,i\geq 1.$$ Such polynomials are orthogonal with respect to the $L^2$-product $$\langle F,G \rangle=\int_0^1 F(s) G(s) w_{a,b}(s)d s,$$ where $w_{a,b}(s)=s^{a-1}(1-s)^{b-1}$ is named the weight function of the basis and is proportional to a beta density in the case of Jacobi polynomials. Any univariate density $f$ supported on $[0,1]$ can be uniquely decomposed on such a basis and therefore there is a unique sequence of real numbers $(\lambda_i)_{i \geq 0}$ such that $$f(s)=w_{a,b}(s)\sum_{i=0}^\infty \lambda_i G_i(s).$$ From the evaluation of $\int_0^1 f(s)\, G_i(s)\,d s$ it follows that each $\lambda_i$ coincides with a linear combination of the first $i$ moments of $S$, specifically $\lambda_i=\sum_{r=0}^i G_{i,r}\gamma_r$. Then, truncate the representation of $f$ in the Jacobi basis at a given level $N$, providing the approximation $$f_N(s)=w_{a,b}(s)\sum_{i=0}^N \left(\sum_{r=0}^i G_{i,r}\mu_r\right) G_i(s).$$ That polynomial approximation is not necessarily a density as it might fail to be positive or to integrate to 1. In order to overcome this problem, one can consider the density $\pi$ proportional to its positive part, thus defined by $\pi(s)\propto\max(f_N(s),0)$. If sampling from $\pi$ is needed, one can resort to the rejection sampler, see for instance Robert & Casella (2).
There is a companion R package called momentify that provides the approximated distribution given moments, and allows to sample from it, available at this link: https://www.researchgate.net/publication/287608666_momentify_R_package, and discussed at this blog post.
Below are two examples with increasing the number of moments involved. Note that the fit is much better for the unimodal density than the multimodal one.
References
(1) Julyan Arbel, Antonio Lijoi, and Bernardo Nipoti. Full Bayesian inference with hazard mixture models. Computational Statistics & Data Analysis 93 (2016): 359-372. arXiv link, journal link.
(2) Christian Robert, and George Casella. "Monte Carlo Statistical Methods Springer-Verlag." New York (2004).