We have found a fast and efficient method for the interpolation of nonstationary seismic data. The method uses the fast generalized Fourier transform (FGFT) to identify the space-wavenumber evolution of nonstationary spatial signals at each temporal frequency. The nonredundant nature of FGFT renders a big computational advantage to this interpolation method. A least-squares fitting scheme is used next to retrieve the optimal FGFT coefficients representative of the ideal interpolated data. For randomly sampled data on a regular grid, we seek a sparse representation of FGFT coefficients to retrieve the missing samples. In addition, to interpolate the regularly sampled seismic data at a given frequency, we use a mask function derived from the FGFT coefficients of the low frequencies. Synthetic and real data examples can be used to examine the performance of the method.