The vibrational spectrum of SO2 up to 15 000 cm(-1) is calculated using a low-storage filter-diagonalization method based on the Chebyshev propagation. The Hamiltonian in the Radau coordinates is expressed on a direct product of one-dimensional discrete variable representation (DVR) grids. The extended symmetry-adapted discrete variable representation (ESADVR) is implemented to accelerate the calculation of the action of kinetic energy operators, and multiple symmetry-adapted autocorrelation functions are obtained from the propagation of a single wave packet. Approximately 1000 vibrational energy levels are identified and some of them are assigned according to the nodal structure of the eigenstates. Comparison with experimental data indicates reasonably good agreement (<1%). The agreement, however, deteriorates with increasing energy, implicating imperfection in the potential energy surface used in the calculation. Statistical analyses indicate that the system is mostly regular in this energy range. There is some evidence of a normal-to-local mode transition at higher energies.