There is the standard argument, using the definition of the inner product; that $\langle f|A|g\rangle =\langle g|A|f\rangle ^{*}$ for a Hermitian operator $A$, given any wave vectors $|f\rangle,~ |g\rangle$.
Also consider the following:
Consider the infinite dimensional one dimensional position space, with a column vector of values of the wave function at discretized points along the $x$-axis.
In the limit of infinitesimal differences between position values, we have $\frac{\mathrm df}{\mathrm dx}\approx \frac{1}{2h} (f(x+h)-f(x-h))$, where $(x+h)$ and $(x-h)$ are the discrete position values just preceding and succeeding the position value $x$, and $h$ is sufficiently small.
Then we might talk of an infinite dimensional matrix representation of ${\rm d/d}x$, where only the two off diagonal "diagonals" adjacent to the actual diagonal has $1/2h$ and $-\:1/2h$ entries. This matrix is skew symmetric. Everywhere else we have $0$ entries.
If we multiplied this matrix with $'\mathrm i'$, this skew symmetry becomes Hermitian, which makes ${\rm i~ d/d}x$ hermitian.
Edit: As pointed out by tparker in an answer below, we get $\langle x|\partial|x^\prime\rangle =\frac{\partial}{\partial x}\langle x|x^\prime \rangle =\frac{\partial}{\partial x}\delta(x-x^\prime)$.
Since we are dealing with a discrete set of points in space here, we must have the normalization given by the Kronecker delta $\delta_{x,x^\prime}$.
Informally, we then have $\partial(\delta_{x,x^\prime})|_{x~=~(x^\prime-h)}\approx\frac{1}{2h}(\delta_{x^\prime,x^\prime}-\delta_{x^\prime-2h,x^\prime})=\frac{1}{2h}$, and also $\partial(\delta_{x,x^\prime})|_{x~=~(x^\prime+h)}\approx-\frac{1}{2h}, ~~\partial(\delta_{x,x^\prime})|_{x~=~x^\prime}\approx 0$, which again gives us a skew-symmetric matrix of the form obtained before. Here we scanned the matrix vertically in a given column, whereas in the previous calculation, we scanned a fixed row horizontally.