Three-dimensional reverse-time migration of teleseismic receiver functions using the phase-shift-plus-interpolation method