Hagedorn functions are carefully constructed generalizations of Hermite functions to the setting of many-dimensional squeezed and coupled harmonic systems. Wavepackets formed by superpositions of Hagedorn functions have been successfully used to solve the time-dependent Schrödinger equation exactly in harmonic systems and variationally in anharmonic systems. To evaluate typical observables, such as position or kinetic energy, it is sufficient to consider orthonormal Hagedorn functions with a single Gaussian center. Instead, we derive various relations between Hagedorn bases associated with different Gaussians, including their overlaps, which are necessary for evaluating quantities nonlocal in time, such as the time correlation functions needed for computing spectra. First, we use the Bogoliubov transformation to obtain the commutation relations between the ladder operators associated with different Gaussians. Then, instead of using numerical quadrature, we employ these commutation relations to derive exact recurrence relations for the overlap integrals between Hagedorn functions with different Gaussian centers. Finally, we present numerical experiments that demonstrate the accuracy and efficiency of * Author to whom any correspondence should be addressed.