What does Mellin inversion "really mean"? Given a function $f: \mathbb{R}^+ \rightarrow \mathbb{C}$ satisfying suitable conditions (exponential decay at infinity, continuous, and bounded variation) is good enough, its Mellin transform is defined by the function
$$M(f)(s) = \int_0^{\infty} f(y) y^s \frac{dy}{y},$$
and $f(y)$ can be recovered by the Mellin inversion formula:
$$f(y) = \frac{1}{2\pi i} \int_{\sigma - i \infty}^{\sigma + i \infty} y^{-s} M(f)(s) ds.$$
This is a change of variable from the Fourier inversion formula, or the Laplace inversion formula, and can be proved in the same way. This is used all the time in analytic number theory (as well as many other subjects, I understand) -- for example, if $f(y)$ is the characteristic function of $[0, 1]$ then its Mellin transform is $1/s$, and one recovers the fact (Perron's formula) that 
$$\frac{1}{2\pi i} \int_{2 - i \infty}^{2 + i \infty} n^{-s} \frac{ds}{s}$$
is equal to 1 if $0 < n < 1$, and is 0 if $n > 1$. (Note that there are technical issues which I am glossing over; one integrates over any vertical line with $\sigma > 0$, and the integral is equal to $1/2$ if $n = 1$.)
I use these formulas frequently, but... I find myself having to look them up repeatedly, and I'd like to understand them more intuitively. Perron's formula can be proved using Cauchy's residue formula (shift the contour to $- \infty$ or $+ \infty$ depending on whether $n > 1$), but this proof doesn't prove the general Mellin inversion formula.
My question is:

What do the Mellin transform and the inversion formula mean? Morally, why are they true?

For example, why is the Mellin transform an integral over the positive reals, while the inverse transform is an integral over the complex plane?
I found some resources -- Wikipedia; this MO question is closely related, and the first video in particular is nice; and a proof is outlined in Iwaniec and Kowalski -- but I feel that there should be a more intuitive explanation than any I have come up with so far.
 A: Have a look at Zagier's appendix:
http://people.mpim-bonn.mpg.de/zagier/files/tex/MellinTransform/fulltext.pdf
It provides a nice description of the Mellin transform when $f(x)$ is sufficiently smooth at $x=0$, and of rapid decay at infinity.
For example, assume $f(x) = \sum_0^\infty a_n x^n$, in some neighbourhood of the origin, and decays rapidly as $x \to \infty$, then its Mellin transform has meromorphic continuation to all of $\mathbb{C}$ with simple poles of residue $a_n$ at $s=-n$, $n=0,1,2,3,\ldots$. This is nicely explained in Zagier's appendix.
So, rate of decay issues aside, shifting the inverse Mellin transform to the left, i.e. letting $\sigma \to -\infty$, picks up the residues of the integrand at s=-n, i.e. $a_n x^n$, i.e. recovers the Taylor expansion about $x=0$ of $f(x)$.
Of course, it only applies to a limited class of functions $f$, but, in many practical examples, this reasoning gives one explanation of why the Mellin inversion formula is true, without resorting to Fourier inversion. 
A: Another property: The (inverse) Mellin transform interchanges $q$-expansions of modular forms and Dirichlet $L$-series.
A: [Some next-day edits in response to comments] As counterpoint to other viewpoints, one can say that Mellin inversion is "simply" Fourier inversion in other coordinates. Depending on one's temperament, this "other coordinates" thing ranges from irrelevancy to substance... The question about moral imperatives for Fourier inversion is addressed a bit below.
[Added: the exponential map $x\rightarrow e^x$ gives an isomorphism of the topological group of additive reals to multiplicative. Thus, the harmonic analysis on the two is necessarily "the same", even if the formulaic aspects look different. The occasional treatment of values (and derivatives) at $0$ for functions on the positive reals, as in "Laplace transforms", is a relative detail, which certainly has a corresponding discussion for Fourier transforms.]
The specific riff in Perron's identity in analytic number theory amounts to (if one tolerates a change-of-coordinates) guessing/discerning an L^1 function on the line whose Fourier transform is (in what function space?!) the characteristic function of a half-line. 
Since the char fcn of a half-line is not in L^2, and does not go to 0 at infinity, there are bound to be analytical issues... but these are technical, not conceptual. 
[Added: the Fourier transform families $x^{\alpha-1}e^{-x}\cdot \chi_{x>0}$ and $(1+ix)^{-\alpha}$, (up to constants) where $\chi$ is the characteristic function, when translated to multiplicative coordinates, give one family approaching the desired "cut-off" effect of the Perron integral. There are other useful families, as well.]
To my taste, the delicacies/failures/technicalities of not-quite-easily-legal aspects of Fourier transforms are mostly crushed by simple ideas about Sobolev spaces and Schwartz' distributions... tho' these do not change the underlying realities. They only relieve us of some of the burden of misguided fussiness of some self-appointed guardians of a misunderstanding of the Cauchy-Weierstrass tradition.
[Added: surely such remarks will strike some readers as inappropriate poesy... but it is easy to be more blunt, if desired. Namely, in various common contexts there is a pointless, disproportionate emphasis on "rigor". Often, elementary analysis is the whipping-boy for this impulse, but also one can see elementary number theory made senselessly difficult in a similar fashion. Supposedly, the audience is being made aware of a "need/imperative" for care about delicate details. However, in practice, one can find oneself in the role of the Dilbertian "Mordac the Preventer (of information services)" [see wiki] proving things like the intermediate value theorem to calculus students: it is obviously true, first, or else one's meaning of "continuous" or "real numbers" needs adjustment; nevertheless, the traditional story is that this intuition must be delegitimized, and then a highly stylized substitute put in its place. What was the gain? Yes, something foundational, but time has passed, and we have only barely recovered, at some expense, what was obviously true at the outset.
On another hand, Bochner's irritation with "distributions theory" was that it was already clear to him that things worked like this, and he could already answer all the questions about generalized functions... so why be impressed with Schwartz' "mechanizing" it? For me, the answer is that Schwartz arranged a situation so that "any idiot" could use generalized functions, whereas previously it was an "art". Yes, sorta took the fun out of it... but maybe practical needs over-rule preservation of secret-society clubbiness?] 
Why should there be Fourier inversion? (for example...) Well, we can say we want such a thing, because it diagonalizes the operator $d/dx$ on the line (and more complicated things can be said in more complicated situations). 
Among other things, this renders "engineering math" possible... That is, one can understand and justify the almost-too-good-to-be-true ideas that seem "necessary" in applied situations... where I can't help but add "like modern number theory". :)
[Added: being somewhat an auto-didact, I was not aware until relatively late that "proof" was absolutely sacrosanct. To the point of fetishism? In fact, we do seem to collectively value insightful conjecture and not-quite-justifiable heuristics, and interesting unresolved ideas offer more chances for engagement than do settled, ironclad, finished discussions. For that matter, the moments that one intuits "the truth", and then begins looking for reasons, are arguably more memorable, more fun, than the moments at which one has dotted i's and crossed t's in the proof of a not-particularly-interesting lemma whose truth was fairly obvious all along. More ominous is the point that sometimes we can see that something is true and works despite being unable to "justify" it. Heaviside's work is an instance. Transatlantic telegraph worked fine despite...]
In other words: spectral decomposition and synthesis. Who couldn't love it?!
[Added: and what recourse do we have than to hope that reasonable operators are diagonalizable, etc? Serre and Grothendieck (and Weil) knew for years that the Lefschetz fixed-point theorem should have an incarnation that would express zeta functions of varieties in terms of cohomology, before being able to make sense of this. Ngo (Loeser, Clucker, et alteri)'s proof of the fundamental lemma in the number field case via model theoretic transfer from the function field case is not something I'd want to have to "justify" to negativists!]
A: Thanks to everyone who answered! A (CW'ed) summary of some of what I learned:
In the first place, I now cheerfully second Greg Martin's recommendation of Chapter 5.1 of Montgomery and Vaughan. It is a rather "lowbrow", very readable treatment. (doesn't prove Mellin inversion in complete generality)
Also, as Matt Young pointed out, for any complex $s$, the function $t \rightarrow t^s$ is a character on $\mathbb{R}^{\times}$. This is a triviality, but the importance of this fact escaped me the first time. The invariant measure on $\mathbb{R}^{\times}$ is $\frac{dx}{x}$, and so the Fourier transform of a function $f$ defined on this group is exactly
$$\int_{x \in \mathbb{R}^{\times}} f(x) x^s \frac{dx}{x},$$
the Mellin transform. Once this is written down, the rest follows mechanically (from change of variables and Fourier inversion).
Thanks to all!
A: As others have pointed out, the Mellin inversion theorem is just the Fourier inversion theorem in disguise for the particular group ${\mathbb R_+}$ with invariant measure $\frac{dx}{x}$.  The goal of the Fourier transform is to express a general function as a linear combination (i.e. integral) of the characters of the group, so that in this basis the operations of translations and all commuting operations will be diagonalized.  For the ${\mathbb R_+}$, these characters look like $x \mapsto x^{-s}$ (the minus sign because of the normalization you chose in the question), and they are unitary (take values in the circle) for imaginary $s$ -- the operation of multiplying characters is just addition in the $s$ variable, so in the inversion formula you have the measure $ds$.  There's also this funny thing about how there are $s$ with positive real part -- this is because in the "physical space" ${\mathbb R_+}$ you're always talking about distributions which are compactly supported away from $0$ when you use this transform.  Let's ignore that.
Since Mellin inversion is a disguised Fourier inversion, the real question is: why is the Fourier inversion formula on ${\mathbb R}$ true?  To me the most convincing answer is the following: we can decompose a general function $f(x) = \int f(y) \delta(x-y) dy$ (this is the definition of $\delta$ but you have to take approximate delta-functions to make this rigorously work like a decomposition), so if we want to express a general function as a combination of the characters $x \mapsto e^{2 \pi i \xi x}$, it suffices to consider the $\delta$ function
$\delta(x) = \int u(\xi) e^{2 \pi i \xi x} d\xi $
One interpretation of this formal idea is that the distributions $\delta(x-y)$ are just like your usual standard basis functions.
Now, observe that because $\delta(x)$ is invariant under multiplication by $e^{2 \pi i \eta x}$ for any $\eta$, the distribution $u(\xi)$ is translation invariant, and therefore must be constant.  After you find the constant, plugging in $\delta(x) = C \int e^{2 \pi i \xi x} d\xi$ into $f(x) = \int f(y) \delta(x-y) dy$ gives the Fourier inversion formula.  Complete, rigorous proofs all follow more or less these lines, but there are many flavors of how you like to phrase it.  Of course, we can write the whole argument with multiplicative characters as well.
Edit:  The above argument assumes uniqueness of the representation, but one can also remark that if there is even a single function $f(x)$ for which $\int f(x) dx \neq 0$ and which can be realized as a linear combination $\int \hat{f}(\xi) e^{2 \pi i \xi \cdot x} d\xi$, then by rescaling, renormalizing and taking a limit, we obtain $\delta(x) = C \lim_{\epsilon \to 0} \epsilon^{-1} f(x/\epsilon)$, leading formally to the formula $\delta(x) = C \int e^{2 \pi i \xi \cdot x} d\xi$.  One common rigorous execution of this philosophy is performed by taking $f$ to be a Gaussian.
A: Two equations that encapsulate the properties of the Fourier and Mellin transforms:
$$\int^{\infty}_{-\infty}{\exp(2 \pi ifx)\exp(-2 \pi ify)df} = \delta(x-y)$$
$$\frac{1}{2\pi i} \int_{\sigma - i \infty}^{\sigma + i \infty} x^{-s} y^{s} ds= \delta(\ln(x)-\ln(y))= y \delta(x-y).$$
The transformations from one equation to the other are obvious. The delta function results are intuitive and an extrapolation of the discrete case for the orthogonality relationships of the characters of character groups. The transform pairs, Plancherel and convolution theorems, and other relations are easy to derive from these two.
(Note that whereas $e^{sz}$ is an eigenfunction for $d/dz$ and so the Laplace/Fourier transforms are appropriate for devising an operator calculus for $f(d/dz)$, $z^s$ is an eigenfunction of $zd/dz$ and so the Mellin transform is more appropriate for $f(zd/dz)$.)
Ramanujan's Master Formula/Theorem (see Wikipedia, particularly the Hardy refer.) gives a somewhat intuitive perspective on the Mellin transform as providing an "interpolation" of the coefficients of the Taylor series of certain classes of functions, as discussed in the intro of "Ramanujan's Master Theorem ..." by Olafsson and Pasquale. E.g.,
$$\int^{\infty}_{0}f(x)\frac{x^{s-1}}{(s-1)!} dx = g(-s)$$  and
$$\frac{1}{2\pi i} \int_{\sigma - i \infty}^{\sigma + i \infty} \frac{\pi}{\sin(\pi s)} g(-s) \frac{x^{-s}}{(-s)!} ds = \sum_{n=0}^{\infty} g(n) \frac{(-x)^{n}}{n!} = f(x)$$
for the transform pairs
$f(x)=\exp(-x)$  and  $g(-s)= 1$ $(\sigma>0)$  and
$f(x)=\frac{1}{1+x}$  and  $g(-s)= (-s)!$ $(0<\sigma<1$  and  $abs(x)<1)$
$f(x)=\exp(-x^2)$  and  $g(-s)= \cos(\pi\frac{ s}{2})\frac{(-s)!}{(-\frac{s}{2})!} = \frac{1}{2}\frac{(\frac{s}{2}-1)!}{(s-1)!}
$ $(\sigma>0)$.
From a similar perspective, the iconic Euler (Mellin) integral for the gamma function for $Real(s) > 0$
$$ \displaystyle \int_0^\infty \frac{t^{s-1}}{(s-1)!} \; e^{-t\;p} \; dt = p^{-s}$$
provides the scaffolding for understanding and utilizing the interplay among the Mellin transform, its inverse, operator calculus, and interpolation.
A natural interpolation of the derivative as the fractional integroderivative of fractional calculus is obtained by using the Mellin transform to interpolate the op coefficients of the op e.g.f. $\displaystyle e^{tD_x} \;,$ i.e., the shift op, for the integer powers of the derivative:
$$\displaystyle \int_0^\infty \frac{t^{s-1}}{(s-1)!} \; e^{-tD_x} \; dt \; H(x) g(x) = D_x^{-s} H(x) g(x) = \int_0^\infty \frac{t^{s-1}}{(s-1)!} \; e^{-tD_x} \; H(x) g(x)\; dt$$
$$\displaystyle = \int_0^\infty \frac{t^{s-1}}{(s-1)!} \; H(x-t) \; g(x-t) dt \; . $$
Then specifically acting on the power function for $\displaystyle \alpha > -1$
$$ \displaystyle \int_0^\infty \frac{t^{s-1}}{(s-1)!} \; H(x-t) \; (x-t)^\alpha dt = \int_0^x \frac{t^{s-1}}{(s-1)!} \; (x-t)^\alpha \; dt $$
$$\displaystyle = \int_0^x \frac{t^{s-1}}{(s-1)!} \; \sum_{k \ge 0} (-1)^k \; x^{\alpha-k} \frac{\alpha!}{(\alpha-k)} \; \frac{t^k}{k!} \; dt = \frac{1}{(s-1)!} \sum_{k \ge 0} (-1)^k \; x^{\alpha-k} \binom{\alpha}{k} \; \frac{t^{s+k}}{s+k} \; |_{t=0}^{x}$$
$$\displaystyle = x^{\alpha + s} \; (-s)! \; \sum_{k \ge 0} \; \binom{\alpha}{k} \; \frac{sin(\pi (s+k))}{\pi (s+k)} = x^{\alpha +s} \frac{\alpha!}{(\alpha+s)!} \; = D_x^{-s} x^\alpha \; .$$
The last summation converges with no restriction on $s$. So, we see that the Mellin transform does indeed interpolate the coefficients of the e.g.f. generated by the binomial theorem expansion $\displaystyle x^{\alpha-k} \frac{\alpha!}{(\alpha-k)}$ to $\displaystyle x^{\alpha+s} \frac{\alpha!}{(\alpha+s)}$ to give an interpolation of the coefficients of the shift op $ D_x^k$ to $ D_x^{-s}$ consistent with fractional calculus.
The same method can be used to interpolate
$$\displaystyle (x \; D_x \;x)^n = x^n D_x^n x^n = x^n \; n!\; L_n(-:xD_x:) , $$
where $ L$ denotes the Laguerre polynomials and $(:xD_x:)^k = x^kD_x^k$ by definition, leading to
$$ \int_0^\infty \frac{t^{s-1}}{(s-1)!} \; e^{-txD_xx} \; dt \; H(x) x^\alpha = (xD_xx)^{-s}\; x^\alpha = \int_0^\infty \frac{t^{s-1}}{(s-1)!} \; \frac{x^\alpha}{(1+xt)^{\alpha+1}} \; dt = x^{\alpha-s} \frac{(\alpha-s)!}{\alpha!} =  x^{-s} D_x^{-s} x^{-s} \; x^\alpha $$
for $ 0 < Real(s) < \alpha +1 \; .$
Or, give the analytic continuation for a Mellin transform related to a class of differential operators encompassing the Witt Lie algebra:
$$ (x^{1+y}D_x)^{-s} \; x^\alpha = \int_0^\infty \frac{t^{s-1}}{(s-1)!} \; H[\frac{x}{(1+y\;t\;x^y)^{1/y}}] \frac{x^\alpha}{(1+y\;t\;x^y)^{\alpha/y}} \; dt $$
$$= H(y) \; x^{\alpha-sy} y^{-s} \frac{(-s+\alpha/y-1)!}{(\alpha/y-1)!} \;+ \; H(-y) \; x^{\alpha+s|y|} |y|^{-s} \frac{(\alpha/|y|)!}{(\alpha/|y|+s)!} \;.$$
A simple way to derive the formulas in your question is by looking at the inverse Mellin transform rep of the Dirac delta function. See my short note on the Inverse Mellin Transform and the Dirac Delta Function. See also some applications in Dirac's Delta Function and Riemann's Jump Function J(x) for the Primes and The Inverse Mellin Transform, Bell Polynomials, a Generalized Dobinski Relation, and the Confluent Hypergeometric Functions.
Edwards in Riemann's Zeta Function in Ch. 10 Fourier Analysis Sec. 10.1 Invariant Operators on R+ and Their Transforms gives a nice, more group-theoretic intro to the Mellin transform in line with other comments in this stream.
