So Cris Luengo, here is a short MATLAB program demonstrating this. (The purpose of the program was to demonstrate something about a discontinuity that causes the interpolated bandlimited function to exceed the samples by several dB. But you can change the input sequence, $x[n]$, if you want.)
Now, the length of the interpolation filter is $N$ and the number of fractional times between each input sample is $R$. This uses a simple Kaiser-windowed $\operatorname{sinc}(\cdot)$ function, but you can design your reconstruction filter however you think is better.
It's a ridiculously long interpolation filter; it looks 2048 samples in the past and 2048 samples into the future to interpolate between any two adjacent samples. So the brickwall is pretty damn bricky. And there are 256 little fractional delays in between each input sample. Both numbers can be increased to as high as you need to make the precision of this reconstruction to be as good as you need.
And the continuous-time $t$ need not be an integer multiple of $\frac{1}{R}$. You can set $R$ very large and then linearly interpolate between adjacent micro-samples and have $t$ have any precision you want. If you don't like linear interpolation, you can use a higher order polynomial like a cubic Hermite spline and get it even more smooth. Ain't smooth enough? Then increase $R$. Brickwall filter not bricky enough? Then increase $N$.
N = 4096; % number of samples used in interpolation calculation
R = 256; % number of equally-spaced fractional sample delays
stop_band_attenuation = 90;
h = zeros(1, N*R); % zero the whole thang
t = linspace(-N/2+1/R, N/2-1/R, N*R-1);
h = [0 sinc(t)/R .* kaiser(N*R-1, 0.1102*(stop_band_attenuation-8.7))'] ; % kaiser windowed sinc()
clear t;
%
% this defines a simple triangle (linear) interpolation function (for testing purposes)
%
% h((N/2-1)*R+1:(N/2)*R+1) = linspace(0.0, 1.0/R, R+1);
% h((N/2)*R+1:(N/2+1)*R+1) = linspace(1.0/R, 0.0, R+1);
%
xLen = 16384;
t_x = linspace(-xLen/2, xLen/2-1, xLen);
x = [repmat([-1 1], 1, xLen/4) repmat([1 -1], 1, xLen/4)];
yLen = R*(xLen+N);
t_y = linspace(-(xLen+N)/2, (xLen+N)/2-1/R, yLen);
y = zeros(1, yLen);
for n = 0:xLen-1
y(R*n+1:R*n+R*N) = y(R*n+1:R*n+R*N) + x(n+1)*R*h;
end
figure(1)
plot(t_x, x, 'ro')
hold on
plot(t_y, y, 'b')
hold off
figure(2)
plot(R*h); % plot impulse response
pause;
freq = linspace(0, R-R/(8*N*R), 8*N*R); % set of frequencies from DC to just under Nyquist
H = fft([h([N*R/2:N*R-1]) zeros(1, 15*N*R+1) h([1:N*R/2-1])]); % fft after zero-padding to extend length by factor of 8
plot(freq, 20*log10(abs(H(1:8*N*R))+1e-16)); % dB by linear freq plot
pause;
semilogx(freq(2:8*N*R), 20*log10(abs(H(2:8*N*R))+1e-16)); % dB by log freq plot, skip DC
pause;
So here is what the brickwall filter looks like. Is it sufficiently bricky for you?
