A data-driven method combined with the formulations of boundary integral equations is developed for the frequency-domain analysis of parametrized acoustic systems, arising from the spatial discretization of the linear Helmholtz equation. The method derives surrogate models for the approximation of frequency response functions at selected field points via the construction of neural networks with radial basis function interpolation. This requires the offline collection of a training database consisting of discrete frequency-parameter samples and their associated response values, whereas the information about system matrices and interference with boundary element kernels are no longer necessary, leading to a matrix-free and non-intrusive nature. To ensure accurate fitting results, a tensor decomposition technique is employed to decouple the naturally formed tensor grid between the frequency and parameter inputs in the database. After establishing two sets of radial basis networks separately, the final reduced-order solutions as a continuous function of the frequency and parameters are represented by the linear combinations of tensor products. In the online phase, the original transfer function at any testing frequency-parameter location can be quickly predicted as direct output from the compact surrogate model. With the analytical solution, the traditional boundary element method and the fast multipole boundary element method as solvers for the generation of training data, two benchmark problems and one complex system are investigated. In all cases, the simplicity, versatility and efficiency of the proposed method are demonstrated through examples of multi-frequency solutions of acoustic systems with parametrized material and geometric properties.