Suppose you have two states in a finite orthonormal (i.e. orthogonal and normalized) basis, $|a\rangle = \sum_nc_n|x_n\rangle$ and $|b\rangle = \sum_n k_n|x_n\rangle$, where $n$ ranges over some finite set. Computing the inner product of these states is straightforward: it's just $\langle a|b\rangle=\sum_n\sum_m c_n^*k_m\langle x_n|x_m\rangle$. Since the basis is orthonormal, we know that $\langle x_n|x_m\rangle=\delta_{mn}$, where $\delta_{mn}$ is the Kronecker delta (https://en.wikipedia.org/wiki/Kronecker_delta), so we can simplify:
$$\langle a|b\rangle = \sum_m\sum_n c_n^*k_m\delta_{mn}=\sum_n c_n^*k_n$$
Now, suppose we want to construct a similar kind of basis for position space (in one dimension). We'll start by dividing up the $x$ axis into a countable number of intervals of length $\Delta x$, choosing as a set of endpoints $x_n=n\Delta x$, where $n\in\mathbb{Z}$. Once we have a set of intervals, we can construct a set of states that are uniformly distributed over each interval (very much like the rectangles usually employed in Riemann summation). Since each such state occupies the interval $[x_m,x_m+\Delta x)$ for some $m$, we can label them by that interval, constructing the set of states $|[x_m,x_m+\Delta x)\rangle$.
To make this set of states an orthonormal basis, we require two things:
- Each state must be normalized: in other words, $\langle [x_m,x_m+\Delta x)|[x_m,x_m+\Delta x)\rangle=1$ for any $m$. This means that the total probability of finding each state in its respective interval must always be 1.
- The states must be orthogonal: in other words, $\langle [x_m,x_m+\Delta x)|[x_n,x_n+\Delta x)\rangle=0$ for any $m\neq n$. This means that the probability of finding each state anywhere else besides its interval must be 0.
We can now specify a position state in this basis as
$$|f\rangle = \sum_n f_n |[x_n,x_n+\Delta x)\rangle$$
For convenience, instead of labeling the coefficients $f_n$, let's instead make them the outputs of a discrete-domain function $f:\{x_n\}_{n\in\mathbb{Z}}\to\mathbb{C}$, so that the state is now written
$$|f\rangle = \sum_n f(x_n)|[x_n,x_n+\Delta x)\rangle$$
Of course, the space of position states that can be constructed using this basis is quite limited, since the probability of finding the state within any of the intervals of size $\Delta x$ that we defined earlier must be constant for any function constructed using this basis. The smaller we make $\Delta x$, the more position states are accessible. In the end, though, what we really want is to construct any arbitrary position state, no matter how much its amplitude varies as a function of position. What this means is that we need to take the limit $\Delta x\to 0$.
With a finite $\Delta x$, over any finite length $L$, there are at most $\frac{L}{\Delta x}$ basis states; in contrast, with a $\Delta x$ approaching zero, there are now an infinite number of basis states over any finite interval. In taking $\Delta x$ to zero, we have turned our countable basis into an uncountable one, and our discrete basis into a continuous one.
Notice that the last expression we constructed looks nearly exactly like the Riemann sum usually used in integration: we're adding together a bunch of "rectangles" with "height" $f(x_n)$ and "width" $|[x_n,x_n+\Delta x)\rangle$. This means that, when we transfer from a discrete to a continuous basis, we can use the notation and machinery that is usually used when working with Riemann sums: namely, the interval $[x_n,x_n+\Delta x)$ turns into $[x_n,x_n+dx)$ for infinitesimal $dx$, and the sum turns into an integral:
$$|f\rangle = \int f(x)|[x,x+dx)\rangle $$
The state $|[x,x+dx)\rangle$ is uniformly distributed over an interval of infinitesimal length (i.e. length approaching zero). This means it's equivalent to a state which is localized entirely at the point $x$, and so we identify it as simply $|x\rangle$. From this, we get the usual expression:
$$|f\rangle = \int f(x)|x\rangle dx$$
If you want to take the inner product of two such states, remember that $\langle x|y\rangle$ is the zero-width limit of the interval states that we defined earlier, and you'll find, by taking a similar limit, that $\langle x|y\rangle = \delta(x-y)$, where $\delta(x-y)$ is the Dirac delta function. This should make the procedure from above directly applicable.
