The traditional dual-mode ultrasound technology has greater potential for wearable implementation compared to the other emerging technologies. During pulse Doppler blood flow velocity estimation, dual-mode ultrasound needs to emit B-mode pulses simultaneously for imaging localization, which requires that B-mode pulses and Doppler pulses share sampling time. This issue can be addressed by using a sparse interval emission method based on single-frequency Doppler pulses. However, common sparse emission arrangements (such as nested emission, coprime emission, etc.) have the disadvantage of a long sampling time window. Especially, there is the problem of insufficient temporal resolution with significant changes in blood flow velocity. In addition, a longer time window contains more blood flow velocity components, which is prone to artifacts caused by sparse sampling, thereby affecting the accuracy of blood flow velocity estimation. Therefore, this paper proposed a novel blood flow velocity estimation method based on multi-frequency pulse sampling. Firstly, a multi-frequency pulse sampling echo model was constructed and derived. It is proved that under the assumption of stationary frequency band attenuation, this mathematical model is equivalent to the spare interval emission method of single-frequency Doppler pulses. That is, the multi-frequency pulse sampling can achieve performance similar to that of the long time window of the single-frequency pulse sampling mode through a shorter time window, thereby improving the temporal resolution and estimation accuracy of the blood flow velocity spectrum. Subsequently, this paper proposed two methods for constructing the blood flow velocity spectrum for multi-frequency pulse sampling mode, namely the BMUSIC algorithm with lower complexity and the VMUSIC algorithm with better artifact suppression effect. The experimental results based on Field Ⅱ simulation data and in vivo data show that, compared with the sparse emission method of single-frequency Doppler pulses, the proposed method in this paper can not only use a shorter time window to improve the temporal resolution of the blood flow velocity spectrum, but also obtain continuous, clear, high-precision, and better artifact suppression effect of blood flow velocity estimation results. In the experiments with in vivo data, due to the limitation of experimental constraints, the proposed method cannot achieve non-integer multiples of frequency, and the performance advantage of the VMUSIC algorithm is not fully demonstrated.