The fast Fourier transform (FFT) based matrix-free ansatz interpolatory approximations of periodic functions are fundamental for efficient realization in several applications. In this work we design, analyze, and implement similar constructive interpolatory approximations of spherical functions, using samples of the unknown functions at the poles and at the uniform spherical-polar grid locations ( jπ N , kπ N ), for j = 1, . . . , N − 1, k = 0, . . . , 2N − 1. The spherical matrix-free interpolation operator range space consists of a selective subspace of two dimensional trigonometric polynomials which are rich enough to contain all spherical polynomials of degree less than N . Using the O(N 2 ) data, the spherical interpolatory approximation is efficiently constructed by applying the FFT techniques (in both azimuthal and latitudinal variables) with only O(N 2 log N ) complexity. We describe the construction details using the FFT operators and provide complete convergence analysis of the interpolatory approximation in the Sobolev space framework that are well suited for quantification of various computer models.We prove that the rate of spectrally accurate convergence of the interpolatory approximations in Sobolev norms (of order zero and one) are similar (up to a log term) to that of the best approximation in the finite dimensional ansatz space. Efficient interpolatory quadratures on the sphere are important for several applications including radiation transport and wave propagation computer models. We use our matrix-free interpolatory approximations to construct robust FFT-based quadrature rules for a wide class of non-, mildly-, and stronglyoscillatory integrands on the sphere. We provide numerical experiments to demonstrate fast evaluation of the algorithm and various theoretical results presented in the article.