This is more of an explanation about why non-isotropic stationary processes aren't often (or ever) really dealt with explicitly than a positive answer.
For a finite Borel measure $F$ on $\mathbb{R}^m$, write (loosely)
$F(k,dk) := F((k_1,k_1+dk_1) \times \dots \times (k_m,k_m+dk_m))$.
Bochner’s theorem essentially states that $B$ is a wide-sense stationary (WSS) covariance function iff
$B(x) = \int_{\mathbb{R}^m} F_B(k,dk) \exp(2\pi i \langle k, x \rangle)$
for some finite Borel measure $F_B$ called the spectral measure.
In one dimension a mild (absolute continuity) restriction and the Radon-Nikodym theorem allow us to write
$f_B(k) := F_B((-\infty, k)) \Rightarrow B(x) =\int_\mathbb{R} F_B(k,dk) \exp(2 \pi i k x) = \int_\mathbb{R} dk f'_B(k) \exp(2 \pi i k x)$,
which implies that $\Rightarrow \hat B = f'_B$, where the Fourier transform is indicated.
In particular, given the (Fourier transform of the WSS) covariance, the spectral measure is also in hand. In higher dimensions, however, the situation is more difficult. Although we can still write $B(x)$ along the lines above, actual computations are hard. To drive this point home, notice that
$d\nu_g(k) = dx \lvert \det \nabla_k g \rvert \overset{g \in C^1(\mathbb{R}^m,\mathbb{R}^m)}{\Leftrightarrow} \nu_g = \lambda \circ g$,
(where $\lambda$ is Lebesgue measure) does not provide a solution for the inverse problem of obtaining the diffeomorphism $g$ from $\nu_g$.
Isotropy makes things much easier, as you've noticed.