Abstract: An interpolation algorithm based on the Krylov subspace projection is proposed to resolve the computational difficulty in multi-frequency acoustic radiation analysis by the boundary element method.It recasts acoustic quantities to be solved into transfer functions of a linear system.Then a matrix-free model reduction formulation is employed to reduce the original system in the frequency band of interest.Interpolation of the acoustic quantities is realized by the transfer functions of the reduced order system.Equally-spacing and adaptive interpolation schemes are presented to calculate the acoustic radiation pressure and power for three model examples.Numerical results demonstrate that the proposed algorithm can significantly speed up the multi-frequency analysis of acoustic radiation problems with a guaranteed computational accuracy.