We study quasinormal modes of black holes in Lovelock gravity. We formulate the WKB method adapted to Lovelock gravity for the calculation of quasinormal frequencies (QNFs). As a demonstration, we calculate various QNFs of Lovelock black holes in seven and eight dimensions. We find that the QNFs show remarkable features depending on the coefficients of the Lovelock terms, the species of perturbations, and spacetime dimensions. In the case of the scalar field, when we increase the coefficient of the third order Lovelock term, the real part of QNFs increases, but the decay rate becomes small irrespective of the mass of the black hole. For small black holes, the decay rate ceases to depend on the Gauss-Bonnet term. In the case of tensor type perturbations of the metric field, the tendency of the real part of QNFs is opposite to that of the scalar field. The QNFs of vector type perturbations of the metric show no particular behavior. The behavior of QNFs of the scalar type perturbations of the metric field is similar to the vector type. However, available data are rather sparse, which indicates that the WKB method is not applicable to many models for this sector.