This study presents an improved approach to common-conversion point stacking of converted body waves that incorporates scattering kernels, accurate and efficient measurement of stack uncertainties, and an alternative method for estimating free surface seismic velocities. To better separate waveforms into the P and SV components to calculate receiver functions, we developed an alternative method to measure near surface compressional and shear wave velocities from particle motions. To more accurately reflect converted phase scattering kernels in the common-conversion point stack, we defined new weighting functions to project receiver function amplitudes only to locations where sensitivities to horizontal discontinuities are high. To better quantify stack uncertainties, we derived an expression for the standard deviation of the stack amplitude that is more efficient than bootstrapping and can be used for any problem requiring the standard deviation of a weighted average. We tested these improved methods on Sp phase data from the Anatolian region, using multiple bandpass filters to image velocity gradients of varying depth extents. Common conversion point stacks of 23,787 Sp receiver functions demonstrate that the new weighting functions produce clearer and more continuous mantle phases, compared to previous approaches. The stacks reveal a positive velocity gradient at 80-150 km depth that is consistent with the base of an asthenospheric low velocity layer. This feature is particularly strong in stacks of longer period data, indicating it represents a gradual velocity gradient. At shorter periods, a lithosphere-asthenosphere boundary phase is observed at 60-90 km depth, marking the top of the low velocity layer.