This paper presents a method for calculating the evolutionary power spectral density (EPSD) of the seismic response of bridges using the convolution summation. With zero initial values, a formula for the dynamic component of the response of bridges to spatially varying seismic ground motions is derived as the convolution summation, by assuming the seismic acceleration to vary linearly between two adjacent time stations. The convolution summation is used for calculating the convolution integral of the dynamic component response factor in the EPSD, in which the constant coefficients are independent of the harmonically modulated excitation. The constant coefficients are obtained by the time-history analysis using triangular unit impulse acceleration excitations. The computational cost of the EPSD depends mainly on the amount of degrees-of-freedom (DOFs) numbers of bridge supports in contact. The corresponding computational scheme is proposed, and its validity is indirectly verified with a single DOF system, by comparing the results obtained with those of the existing methods. Finally, a three-span continuous rigid-frame bridge is taken as a case study to illustrate the applicability and effectiveness of the proposed scheme.