The coefficient matrices of conventional boundary element method (CBEM) are dense and fully populated. Special techniques such as hierarchical matrices (H-matrices) format are required to extent its ability of handling large-scale problems. Adaptive cross approximation (ACA) algorithm is a widely adopted algorithm to obtain the H-matrices. However, the accuracy of the ACA boundary element method (ACABEM) cannot be adjusted by changing the tolerance [Formula: see text] when it exceeds a certain value. In this paper, the degenerate kernel approximation idea for the low-rank matrices is developed to build a fast BEM for acoustic problems by exploring the multipole expansion of the kernel, which is referred as the multipole expansion H-matrices boundary element method (ME-H-BEM). The newly developed algorithm compresses the far-field submatrices into low rank submatrices with the expansion terms of Green’s function. The obtained H-matrices are applied in conjunction with the generalized minimal residual method (GMRES) to solve acoustic problems. Numerical examples are carefully set up to compare the accuracy, efficiency as well as memory consumption of the CBEM, ACABEM, fast multipole boundary element method (FMBEM) and ME-H-BEM. The results of a pulsating sphere indicate that the ME-H-BEM keeps both storage and operation logarithmic-linear complexity of the H-matrices format as the ACABEM does. Moreover, the ME-H-BEM can achieve better convergence and higher accuracy than the ACABEM. For the analyzed complicated large-scale model, the ME-H-BEM with appropriate number of expansion terms has an advantage in terms of efficiency as compared with the ACABEM. Compared with the FMBEM, the ME-H-BEM is easier to be implemented.