The time-dependent numerical renormalization group (TDNRG) method [Anders et al. Phys. Rev. Lett. 95, 196801 (2005)] offers the prospect of investigating in a non-perturbative manner the time-dependence of local observables of interacting quantum impurity models at all time scales following a quantum quench. Here, we present a generalization of this method to arbitrary finite temperature by making use of the full density matrix approach [Weichselbaum et al. Phys. Rev. Lett. 99, 076402 (2007)]. We show that all terms in the projected full density matrix ρ i→ f = ρappearing in the time-evolution of a local observable may be evaluated in closed form at finite temperature, with ρ +− = ρ −+ = 0. The expression for ρ −− is shown to be finite at finite temperature, becoming negligible only in the limit of vanishing temperatures. We prove that this approach recovers the short-time limit for the expectation value of a local observable exactly at arbitrary temperatures. In contrast, the corresponding long-time limit is recovered exactly only for a continuous bath, i.e. when the logarithmic discretization parameter Λ → 1 + . Since the numerical renormalization group approach breaks down in this limit, and calculations have to be carried out at Λ > 1, the long-time behavior following an arbitrary quantum quench has a finite error, which poses an obstacle for the method, e.g., in its application to the scattering states numerical renormalization group method for describing steady state non-equilibrium transport through correlated impurities [Anders, Phys. Rev. Lett. 101, 066804, (2008)]. We suggest a way to overcome this problem by noting that the time-dependence, in general, and the long-time limit, in particular, becomes increasingly more accurate on reducing the size of the quantum quench. This suggests an improved generalized TDNRG approach in which the system is time-evolved between the initial and final states via a sequence of small quantum quenches within a finite time interval instead of by a single large and instantaneous quantum quench. The formalism for this is provided, thus generalizing the TDNRG method to multiple quantum quenches, periodic switching, and general pulses. This formalism, like our finite temperature generalization of the single-quench case, rests on no other approximation than the NRG approximation. The results are illustrated numerically by application to the Anderson impurity model.