Given samples of an everywhere non-negative or positive-valued continuous-time signal band-limited to half the sampling frequency, is there some practically applicable way to interpolate it so that there is no danger of producing negative values? The approach should be scalable to arbitrary accuracy so basic linear interpolation won't do. By arbitrary accuracy I mean an arbitrary small difference between the band-limited continuous-time signal and the interpolation.
Perfect sinc interpolation would work because negative contributions would get cancelled by positive contributions. Ideal sinc interpolation is not possible in practice so other filters are used. This could give negative values because usually the filter impulse response has negative values. For example cubic Hermite spline piece-wise polynomial interpolation gives this impulse response:
Figure 1. Cubic Hermite spline impulse response has negative values.
Lagrange polynomials used as splines would also have negative values in the impulse response. B-splines have everywhere non-negative impulse responses, but their pass-band frequency responses deviate increasingly from the brick-wall frequency response of sinc interpolation with increasing B-spline order. The B-spline impulse response tends to a Gaussian function as the order is increased. A Gaussian function in time domain is a Gaussian function in frequency domain, so not a perfect brick-wall filter.
The obvious way to fix the problem would be to 1) take the absolute value of the interpolation output or 2) set all negative output values to zero, but I'm wondering if there is some better (inherently non-negative) method, because both 1) and 2) would create discontinuities in the derivative of the signal. This also creates spectral images.
I don't know how bad the problem stated in this question really is, because if the band-limited signal is non-negative, this also rules out some contradictory discrete signals like $[\dots, 0, 0, 0, 1, 0, 0, 0, \dots].$ The corresponding non-negative continuous signal would then have to be sinc, which is partially negative, which is the contradiction.
Shay Maymon has investigated "bandlimited square-roots" in his 2011 thesis Sampling and Quantization for Optimal Reconstruction. The idea there is to find (not easy) another band-limited signal (complex and at half the bandwidth) the squared magnitude of which equals the wanted band-limited signal. Interpolating that other signal and multiplying the result by its complex conjugate is never negative. That reminds me of this question about interpolation of the magnitude of discrete Fourier transform (DFT).