Getting the kernel of the operator
There are lots of kinds of operators. Linear, non-linear. Some linear operators can be written with a kernel $k(x,y)$, for example
$$ (A_kf)(x):=\int k(x,y)f(y)dy $$
where the integral can go over different areas, curves, or whatever. Since you are talking about convolution operators, I assume that the operators can be written as a convolution
$$ (A_kf)(x):=(k*f)(x):=\int_\Omega k(x-y)f(y)\,dy\,, $$
where $\Omega$ is the set of points where $k$ is not zero. In the discrete case, the integral becomes a summation:
$$ (A_kf)(x)=\sum_{y\in\Omega} k(x-y)f(y). $$
If you take the function
$$ \delta_{x_0}(x):=\begin{cases}1 & \text{for $x=x_0$,} \\ 0 & \text{otherwise} \end{cases} $$
for some arbitrary pixel position $x_0$ then you can get the kernel of the convolution operator by
$$ (A_k\delta_{x_0})(x_0+x)
= \sum_{y\in\Omega} k(x+x_0-y)\underbrace{\delta_{x_0}(y)}_{=\begin{cases}\text{$1$, if $y=x_0$}\\ \text{$0$, otherwise}\end{cases}}
= k(x) $$
as you can easily check. So in image-processing take a black image (filled with zeroes) and set one pixel to $1$. Then apply the operator and you'll get the operator kernel by the above formula. Sometimes there are boundary effects, so make sure the pixel position $x_0$ is far enough away from the boundary of the image.
Getting the Fourier transform of the kernel
As soon as you have your kernel you can throw a 2d discrete Fourier transformation algorithm at the image data and you'll get it.
There's another way to look at all this: A convolution in space coincides with a point-wise multiplication in the frequency domain. To get the kernel in the Fourier domain, i.e. in the frequency space, you can take an image that is $1$ everywhere in the Fourier space. Then apply the inverse transform to this image, then the operator and the the forward Fourier transform. The result is the Fourier transformed kernel.
The fun this is, that the inverse Fourrier transform applied to the image with $1$ everywhere yields $\delta_0$ modulo a multiplicative constant. So this way of looking at it gives you the same result as the first method. (This is only true, if the operator does not behave strangely at the border of the image, but wraps around, i.e. behaves as if the images would be periodically continued.)
Best Answer
To take it one step further than DisintegratingByParts's answer:
The Fourier transform $\mathcal F$ is unitary on $L^2({\mathbb R})$ with spectrum $\{1,i,-1,-i\}$. This means that $\mathcal F$ effectively behaves like a finite-dimensional diagonal operator. In particular, (basically) any square root of $\mathcal F$ can be written as a polynomial of degree $3$ in $\mathcal F$ itself. These polynomials can be deduced by Lagrange interpolation and equal
$$p(x) = \pm \frac{1}{4} (x+1)(x^2+1) \pm \frac{1+i}{4 i \sqrt{2}} (x^2-1)(x+i) \pm \frac{i}{4} (x-1)(x^2+1) \pm \frac{-1+i}{4 i \sqrt{2}} (x^2-1)(x-i).$$
So the answers are:
On $\mathbb R$, the basis functions are Hermite functions like DisintegratingByParts mentioned. In higher dimensions, $L^2(\mathbb R^d)$ can be decomposed into invariant subspaces whose angular parts are given by spherical harmonics. The book of Stein and Weiss has a very nice exposition.
Generally convolution, pointwise multiplication, differentiaion, and multiplication by monomials will not transform into anything nice simply because the different terms in $p({\mathcal F})$ will transform differently. (It's easy to check that in none of the cases do coefficients magically vanish to produce unexpected symmetry).
Because it's a polynomial in $\mathcal F$, it rather directly extends to distributions.