Seismic data interpolation is essential in a seismic data processing workflow, recovering data from sparse sampling. Traditional and deep-learning-based methods have been widely used in the seismic data interpolation field and have achieved remarkable results. In this paper, we develop a seismic data interpolation method through the novel application of diffusion probabilistic models (DPMs). DPM transforms the complex end-to-end mapping problem into a progressive denoising problem, enhancing the ability to reconstruct complex situations of missing data, such as large proportions and large-gap missing data. The interpolation process begins with a standard Gaussian distribution and seismic data with missing traces and then removes noise iteratively with a UNet trained for different noise levels. Our DPM-based interpolation method allows interpolation for various missing cases, including regularly missing, irregularly missing, consecutively missing, noisy missing, and different ratios of missing cases. The ability to generalize to different seismic data sets is also discussed in this paper. Numerical results of synthetic and field data indicate satisfactory interpolation performance of the DPM-based interpolation method in comparison with the f- x prediction filtering method, the curvelet transform method, the low-dimensional manifold method (LDMM), and the coordinate attention-based UNet method, particularly in cases with large proportions and large-gap missing data. Diffusion is all we need for seismic data interpolation.