Based on the theoretical description of Position-Position-Velocity (PPV) statistics in Lazarian & Pogosyan (2000), we introduce a new technique called the Velocity Decomposition Algorithm (VDA) in separating the contribution of turbulent velocity from density fluctuations. Using MHD turbulence simulations, we demonstrate its promise in recovering the velocity caustics in various physical conditions and, in conjunction with the Velocity Gradient Technique (VGT), its prospects in accurately tracing the magnetic field based on pure velocity fluctuations. Employing the theoretical framework developed in Lazarian & Pogosyan (2004), we find that for localized clouds, the velocity fluctuations are most prominent at the wing part of the spectral line, and they dominate the density fluctuations. The same velocity dominance applies to extended HI regions with galactic rotations. Our numerical experiment demonstrates that velocity channels arising from the cold phase of atomic hydrogen (HI) are still strongly affected by the velocity caustics in small scales. We apply the VDA to the HI GALFA-DR2 data corresponding to the high-velocity cloud HVC186+19-114 and high latitude galactic diffuse HI data. Our study confirms the crucial role of velocity caustics in forming linear structures observed within PPV cubes. We discuss the implications of the VDA for both magnetic field studies and predicting polarized galactic emission that acts as the foreground for the Cosmic Microwave Background (CMB) studies. Besides, we address the controversy related to the nature of the filaments in HI channel maps and explain the importance of velocity caustics in the formation of structures in PPV data cubes. The VDA method will allow astronomers to obtain velocity caustics from almost every piece of spectroscopic PPV data and allows direct investigation of the turbulent velocity field in observations.