Abstract: Numerical simulation of the large-scale high frequency sound field is a computational challenging task. To solve the difficulty that the traditional boundary collocation methods are not easy to be applied to large-scale problems thanks to the resulting large-scale fully-populated matrix, a dual-level singular boundary method is proposed in this study. By introducing a dual level structure, the fully-populated matrix is transformed to a large-scale locally supported sparse matrix on fine mesh. The bottleneck of excessive storage requirements and a large number of operations encountered by the traditional singular boundary method is hereby avoided. Secondly, the method uses only coarse mesh nodes to evaluate far-field contributions, and it is a kernel-independent algorithm. In comparison with the fast multipole method, the dual-level singular boundary method performs higher adaptability and flexibility. In addition, the dual level structure plays a role of preconditioner, which makes the method is very efficient for solving matrix with large scale, high rank and high condition number. In scattering sphere example, the dual-level singular boundary method simulates well the acoustic scattering problem with up to dimensionless wavenumber of 160 when the number of degrees of freedom is taken as 100 000. In the benchmark human head sound scattering, the dual-level singular boundary method using 80 000 degrees of freedom performs 78.13% faster than the COMSOL, and it is noted that the computational frequency is up to 25 kHz, which is far beyond the limit of hearing of human ear.