  We present a novel technique for the calculation of dynamical correlation
functions of quantum impurity systems in equilibrium with Wilson's numerical
renormalization group. Our formulation is based on a complete basis set of the
Wilson chain. In contrast to all previous methods, it does not suffer from
overcounting of excitation. By construction, it always fulfills sum rules for
spectral functions. Furthermore, it accurately reproduces local thermodynamic
expectation values, such as occupancy and magnetization, obtained directly from
the numerical renormalization group calculations.
