Concerning $\tilde{x}$, it is self-adjoint (the proof is easy) if it is defined on its natural domain $D(\tilde x)$. For $\gamma_n(x):= (n+1/2)\delta$ if $n \delta \leq x <(n+1)\delta$ and $n \in \mathbb Z$
$$D(\tilde{x}) = \left\{\psi\in
L^2(\mathbb R) \:\left|\:\int_{\mathbb R} |\gamma_n(x)\psi(x)|^2 dx < +\infty
\right.\right\}$$
The spectrum is $\sigma(\tilde{x})= \{\delta (n+1/2)\:|\: n \in \mathbb Z\}$ as expected.
The fact that your claimed approximated momentum operator $\tilde{p}$ makes any sense as it stands is questionable.
The point is that $\tilde{p}$ defined this way is not self-adjoint but only symmetric, while observables must be self-adjoint in order to exploit the spectral calculus technology.
To obtain a symmetric operator (i.e. Hermitean and densely defined) I think that a possibility is to define its domaine $D(\tilde{p})$ as a space of $C^1$ functions (that vanish in $x_n= \delta n$ for all $n \in \mathbb Z$ with their first derivative for a reason I discuss shortly) also requiring that these functions rapidly vanish with the first derivative for $|x|\to \infty$.
It should be possible to relax the smoothness condition using weak derivatives, but the situation does not seem to change remarkably.
Consider the anti-linear operator $(C\psi)(x):= \overline{\psi(-x)}$, the bar denoting the complex conjugation. It is norm preserving and $CC=I$, so it is a conjugation. It seems to me that $C\tilde{p}= \tilde{p}C$ with the domain I introduced above.
Therefore, in view of a theorem due to von Neumann, $\tilde{p}$ admits some self-adjoint extension. A careful analysis of the defect indices of $\tilde{p}$ would classify these extensions. Therefore we are not authorized to say that $\tilde{p}$ has the meaning of an observable, we have to choose a self-adjoint extension of it. Unfortunately, the spectrum depends on this choice. It is by no means obvious what the spectrum of a self-adjoint extension of $\tilde{p}$ is.
What I am saying is that since we do not know the spectrum of the claimed approximated momentum we cannot say how (if) this claimed approximation works as soon as $\delta \to 0$.
It is possible to prove that if the domain is chosen as $C_0^\infty(\mathbb R)$ or ${\cal S}(\mathbb R)$ then there is a unique self-adjoint extension and coincides with the standard one. But we cannot make this choice in view of the form of $\tilde{x}$.
In fact, when computing $[\tilde{x}, \tilde{p}]=0$
we are assuming that the domain of $\tilde{p}$ is invariant under the action of $\tilde{x}$. This indeed happens if $D(\tilde{p})$ is made of the $C^1$ functions rapidly vanishing as $|x|\to \infty$ also vanishing at each $\delta n$ with their first derivative. The last condition eliminates the discontinuities introduced by the action of $\tilde{x}$ giving rise to a function in the same domain.
We conclude that $[\tilde{x}, \tilde{p}]=0$ holds on a domain which does not fix a unique self-adjoint extension of $\tilde{p}$ (and its spectrum consequently) as it would be if this domain were $C_0^\infty(\mathbb R)$ or ${\cal S}(\mathbb R)$.
An alternative definition of $\tilde{p}$ is obtained re-defining its domain $D(\tilde{p})$ by including the (vanishing sufficiently fast at infinity) functions $C^1$ in each open interval $(n\delta, (n+1)\delta)$ that are separately periodic in each such interval. The values these functions attain in the set of points $\delta n$ is irrelevant as this set has zero measure and $L^2$ does not care of zero measure set: We can safely assume that $\tilde{p}\psi$ vanishes thereon by definition as already assumed by the OP. This only concerns mathematics while, physically speaking, this assumption has devastating implications.
With this definition $\tilde{p}$ become essentially self-adjoint, i.e. it has a unique self adjoint extension. The spectrum should be $\sigma(\tilde{p}_{ext}) = \{ \frac{2\pi m}{\delta} \hbar\:|\: m \in \mathbb Z\}$, since in each said interval we find the standard momentum operator on a segment with periodic boundary conditions (to be completely sure I should check some conditions but I am reasonably confident)
With this definition you find $[\tilde{x}, \tilde{p}]=0$ on $D(\tilde{p})$.
The problem with this construction is that $\tilde{p}$ does not generate space displacements of the wavefunctions, $(e^{-ix_0\tilde{p}}\psi)(x)= \psi(x- x_0)$. And it does not seem that the situation improves as soon as $\delta \to 0$. Instead $\tilde{p}$ generates "periodic" displacements in each interval $(n\delta, (n+1)\delta)$ (as it were a circle). For $\delta \to 0$ these periodic displacements become denser and denser but they have nothing to do with geometric translations along $x$, as the ones generated by the true momentum operator.
