Waveform decomposition techniques are commonly used to extract attributes of targets from light detection and ranging (LiDAR) waveforms. The conventional Gaussian decomposition (GD) cannot deal with system waveforms (SWs) with non-Gaussian shapes, whereas the recently proposed B-spline-based decomposition method holds an assumption of similarity transformation. We present a multi-Gaussian decomposition (MGD) algorithm that employs a Gaussian mixture model (GMM) to represent the SW. Compared with the GD algorithm, the MGD algorithm exploits the characteristic of the SW using the GMM and hence can fit the received waveforms better than the GD algorithm. In contrast to the B-spline-based method, the proposed algorithm holds a more reasonable assumption that a received waveform is a convolution result of the SW with the target response, which accords with the ranging principle better. The MGD algorithm was validated using the experimental waveforms with negative tails collected by our self-designed LiDAR system. The GD algorithm and the B-spline-based decomposition were also introduced and studied for comparison. The experimental results show that the GMM with six components fit the SW with an acceptable residual. The experimental results also show that the MGD algorithm produced better decomposition results than the other two algorithms in terms of range retrieval, whereas the B-spline-based decomposition showed the best performance with regards to the root-mean-square-error for waveform fitting.