Fault diagnosis of rolling bearings under variable speed conditions is a challenging task since the vibration signal exhibits time-varying non-stationary characteristics and is usually contaminated by strong noise. Most of the current researches employ the adaptive filtering or signal decomposition methods to obtain the impulse signals caused by the bearing fault before feature extraction, which, however, are not capable of removing the in-band noise. To address this issue, a novel method called non-stationary harmonic summation (NHS) is proposed based on the fact that the repetitive impulses caused by the bearing fault consist of a series of equally-spaced harmonics in the frequency domain. Firstly, the harmonic characteristics are theoretically analyzed and the results show that the impulses contain non-stationary harmonics with a time-varying spacing frequency (i.e., the fault characteristic frequency) under variable speed conditions. Next, according to the harmonic characteristics, an efficient algorithm combining the parameterized demodulation with the adaptive chirp mode decomposition is developed to extract the non-stationary harmonics and then summate these harmonics to reconstruct the repetitive impulses for the fault feature extraction. Since the NHS elaborately exploits the intrinsic harmonic structure of the impulse signals, the noise can be fully removed and the reconstructed signal is free of side-band interference caused by complex amplitude modulation. Both simulated and experimental signals are considered to demonstrate the advantages of the NHS for bearing fault diagnosis under variable speed conditions.