Does anyone have a (semi-)intuitive explanation of why momentum is the Fourier transform variable of position?
(By semi-intuitive I mean, I already have intuition on Fourier transform between time/frequency domains in general, but I don't see why momentum would be the Fourier transform variable of position. E.g. I'd expect it to be a derivative instead.)
The other answer are all correct, but I want to give another outlook, which I believe is the most self-consistent one. Clearly, we have to start somewhere. My starting point is to define momentum as the quantity that is conserved as a consequence of translational symmetry.
In quantum mechanics the translation operator for a single spin-less particle acts on the states as $$ (T_a \psi)(x)=\psi(x-a) $$ and its generator is$^1$ $$ G=i\hbar\frac{d}{dx} \quad\iff\quad T_a=e^{iaG/\hbar}. $$
If the system is translationally symmetric the generator $G$ is conserved, as follows from general arguments in quantum mechanics. We conventionally consider its opposite as the momentum operator, that is clearly conserved as well.
Now, the eigenstates of the momentum operator are the states $\lvert p\rangle$ whose wavefunctions are$^2$ $$ \langle x\vert p\rangle=e^{ipx/\hbar} $$ and therefore the probability of a generic state $\vert\psi\rangle$ to have momentum $p$ is $$ \langle p\vert\psi\rangle=\int\!dx\,\langle p\vert x\rangle\langle x\vert\psi\rangle =\int\!dx\,e^{-ipx/\hbar}\psi(x). $$ This result says that the probability density of the momentum is the Fourier transform of the probability density of the position$^3$. This is how position and momentum are related through the Fourier transform.
useless nitpicking:
$^1$a way to see this is to compare the expansions \begin{align} (T_a \psi)(x) &= \psi(x-a) = \psi(x)-a\frac{d\psi}{dx}(x) + o(a^2)\\ &= (e^{iaG/\hbar}\psi)(x) = \psi(x) + \frac{ia}{\hbar}(G\psi)(x) + o(a^2) \end{align}
$^2$this is because in position representation the eigenvalue equation for the momentum reads $$ -i\hbar\frac{d\phi_p}{dx}(x)=p\phi_p(x) $$ that has the (unique) aforementioned solution (note that they are not normalizable so they are actually generalized eigenstates)
$^3$a bit loosely speaking, since the probabilities are their $|\cdot|^2$
Momentum is not the Fourier transform of position.
In the position representation, position is the operator of multiplication by $x$, whereas momentum is a multiple of differentiation with respect to $x$. These observables (operators) are not Fourier transforms of each other.
In the momentum representation, momentum is the operator of multiplication by $p$, whereas position is a multiple of differentiation with respect to $p$. These observables (operators) are not Fourier transforms of each other.
The reason why these representations are appropriate for position and momentum is the fact that in both representations, the commutators satisfy the canonical commutation relations, the quantum analogue of the Poisson bracket relation $\{p,q\}=1$.
The Fourier transform comes in only as the means to switch from the position representation to the momentum representation or conversely. The reason is that apart from a factor, differentiation of the Fourier transform of a function $\psi$ is equivalent to multiplication of $\psi$, and differentiation of $\psi$ is equivalent to multiplication of the Fourier transform of $\psi$.
For me this is intrinsically and quite simply a generalization of de Broglie's relation, $$p=\frac{h}{\lambda}=\hbar k.$$ Of course, it this form it only holds for plane-wave kinda wave/particles. In the general case, as Schrödinger posits, the particle is described by some function $\psi(x)$ of position, which is nonzero in some definite range of space. Because you don't have a plane wave any more, you must have some range of wavelengths/momenta in play; if only you could translate your position-function into a function of wavelength (i.e. of momentum) then you'd be sorted. Luckily, this is what a Fourier transform does.
An intuitive explanation based on the time/freq understanding of the Fourier transform comes about by considering the pair, $x$ and $p/\hbar$ following the de Broglie relation Emilio points out. Indeed $p/\hbar$ has units of $1/distance$. Thus, you can consider $t\rightarrow x$ and $\omega\rightarrow k=p/\hbar$.
From here we can align the variable found in time/freq Fourier theory with position/momentum Fourier theory: $$\begin{array}{lll} \textrm{Time/freq} & \textrm{Spc/momentum} \\ \hline t_n(s)&x_n&[n^{th} \textrm{ sample}]\\ N & M&[\textrm{Number of samples}]\\ \omega_m \left(\frac{rad}{s}\right)& k_m=\frac{p}{\hbar} \left(\frac{rad.}{m}\right)&\textrm{FT ang. freq. sample}\\ f_i=\frac{\omega_i}{2\pi} \left(\frac{1}{s}\right)& \tilde{\nu}_i=\frac{k_i}{2\pi} \left(\frac{1}{m}\right)&\textrm{FT freqency samples}\\ f(x_n)&f(t_n)&\textrm{Given function}\\ F(\omega_m)=\sum_n^{N-1}f(t_n)e^{-i\omega_mt_n} &\begin{array}{ll}F(k_m)=\sum_n^{N-1}f(x_n)e^{-ik_mx_n}\end{array}&\textrm{Fourier transform}\\ \end{array}$$ $k_i$ are called (angular) wave numbers and if you consider ordinary frequency, then $\tilde\nu=k_i/2\pi=p/h$ is the spectroscopic wave number (analogous to $f_i=\omega_i/2\pi$).
Momentum is defined by how a solution of the wave equation evolves dynamically, so the relationship given by the fourier transform is not general, but limited to certain evolution laws. The Schroedinger equation of the free particle has certain solutions of the form $\exp(i(\mathbf{k}.\mathbf{r}-\omega t))$ that do nothing but translate in space with a constant velocity. Their velocity and momentum is proportional to $\mathbf{k}$. And these solutions are also the only solutions that have a definite momentum. All other solutions are linear combinations and therefore momentum superpositions.
If you want to know which momentum contributions are present in a given wave function, you have to expand that function in terms of the solutions of definite momentum. Fortunately, these solutions form an orthogonal basis of the Hilbert space so that the expansion becomes a simple inner product, namely $c(\mathbf{k},t) = \int_{-\infty}^\infty \exp(i(\mathbf{k}.\mathbf{r}-\omega t))^* \psi(\mathbf{r})d\mathbf{r}$. We can now apply the complex conjugation and move the time dependency out of the integral: $c(\mathbf{k},t)=\exp(i\omega t)\int_{-\infty}^{\infty}\exp(-i\mathbf{k}.\mathbf{r})\psi(\mathbf{r})d\mathbf{r}$. We have not assumed any time dependency of $\psi$, and so we can just assume that $\psi$ refers to a state at $t=0$ so that we get rid of the leading factor and have $c(\mathbf{k})=\int_{-\infty}^{\infty}\exp(-i\mathbf{k}.\mathbf{r})\psi(\mathbf{r})d\mathbf{r}$, which is the Fourier transform of $\psi(\mathbf{r})$.
