We investigate a variably distributed-order time-fractional wave partial differential equation, which could accurately model, e.g., the viscoelastic behavior in vibrations in complex surroundings with uncertainties or strong heterogeneity in the data. A standard composite rectangle formula of mesh size σ is firstly used to discretize the variably distributed-order integral and then the L-1 formula of degree of freedom N is applied for the resulting fractional derivatives. Optimal error estimates of the corresponding fully-discrete finite element method are proved based only on the smoothness assumptions of the data. To maintain the accuracy, setting σ = O(N−1) leads to O(N3) operations of evaluating the temporal discretization coefficients. To improve the computational efficiency, we develop a novel time-stepping scheme by expanding the fractional kernel at a fixed fractional order to decouple the fractional operator from the variably distributed-order integral. Only O(logN) terms are needed for the expansion without loss of accuracy, which consequently reduce the computational cost of generating coefficients from O(N3) to O(N2 logN). Optimal-order error estimates of this time-stepping scheme are rigorously proved via novel and different techniques from the standard analysis procedure of the L-1 methods. Numerical experiments are presented to substantiate the theoretical results.