In this paper, we propose a parallel-in-time algorithm for approximately solving parabolic equations. In particular, we apply the $k$-step backward differentiation formula, and then develop an iterative solver by using the waveform relaxation technique. Each resulting iteration represents a periodic-like system, which could be further solved in parallel by using the diagonalization technique. The convergence of the waveform relaxation iteration is theoretically examined by using the generating function method. The approach we established in this paper extends the existing argument of single-step methods in Gander and Wu [Numer. Math., 143 (2019), pp. 489--527] to general BDF methods up to order six. The argument could be further applied to the time-fractional subdiffusion equation, whose discretization shares common properties of the standard BDF methods, because of the nonlocality of the fractional differential operator. Illustrative numerical results are presented to complement the theoretical analysis.