Summary: We consider the problem of the reconstruction of a continuous-time function \(f(x)\in{\mathcal H}\) from the samples of the responses of \(m\) linear shift-invariant systems sampled at \(1/m\) the reconstruction rate. We extend Papoulis' generalized sampling theory in two important respects. First, our class of admissible input signals (typ. \({\mathcal H}=L_2\)) is considerably larger than the subspace of band-limited functions. Second, we use a more general specification of the reconstruction subspace \(V(\varphi)\), so that the output of the system can take the form of a band-limited function, a spline, or a wavelet expansion. Since we have enlarged the class of admissible input functions, we have to give up Shannon and Papoulis' principle of an exact reconstruction. Instead, we seek an approximation \(\widetilde f\in V(\varphi)\) that is consistent in the sense that it produces exactly the same measurements as the input of the system. This leads to a generalization of Papoulis' sampling theorem and a practical reconstruction algorithm that takes the form of a multivariate filter. In particular, we show that the corresponding system acts as a projector from \({\mathcal H}\) onto \(V(\varphi)\). We then propose two complementary polyphase and modulation domain interpretations of our solution. The polyphase representation leads to a simple understanding of our reconstruction algorithm in terms of a perfect reconstruction filter bank. The modulation analysis, on the other hand, is useful in providing the connection with Papoulis' earlier results for the band-limited case. Finally, we illustrate the general applicability of our theory by presenting new examples of interlaced and derivative sampling using splines.