Sinc interpolation in spatial domain
The sinc function actually represents an ideal (brickwall) lowpass filter that's used to complete the interpolation process after the data has been expanded (zero stuffed) properly. So let me outline the time domain approach here:
Assume you have data samples $x[n]$ of length $N$, and you want to upsample this by the integer factor of $L$, yielding a new interpolated data $y[n]$ of length $M = L \times N $ samples.
The first stage is to zero stuff the input $x[n]$; i.e., expand it by $L$ by the expander block :
$$ x[n] \longrightarrow \boxed{ \uparrow L } \longrightarrow x_e[n] $$
where $x_e[n]$ is related to $x[n]$ by the following:
$$ x_e[n] = \begin{cases} { x[n/L] ~~~,~~~ n = 0, \pm L, \pm 2L... \\ ~~~ 0 ~~ ~~~,~~~\text{otherwise} } \end{cases} $$
Then, to complete the interpolation process and fill in the empty (zeroed) samples of $x_e[n]$, one has to lowpass filter $x_e[n]$ by an ideal lowpass filter $h[n]$ with the following frequency domain definition $H(\omega)$ :
$$ H(\omega) = \begin{cases} { ~ L ~ ~~~,~~~ |\omega| &lt; \frac{\pi}{L}  \\ ~ 0 ~ ~~~,~~~\text{otherwise} } \end{cases} $$
The impulse response $h[n]$ of this ideal filter is computed by the inverse discrete-time Fourier transform of $H(\omega)$ and is given by
$$ h[n] = L \frac{ \sin( \frac{\pi}{L} n) } {\pi n} $$
This is an infitely long and non-causal filter, and thus cannot be implemented in this form. (See Hilmar's comments) Practically it's truncated and weighted by a window function, for example by a Hamming or Kaiser window.
The following MATLAB / OCTAVE code represents designing the filter and applying it into data in time domain:
L = 5;      % interpolation factor
N = 500;    % data length
 
x = hamming(N)'.*randn(1,N);  % generate bandlimited data...

% expanded signal
xe = zeros(1,N*L);
xe(1:L:end) = x;          % generate th expanded signal 

% interpolation sinc filter
n = -32:32;               % timing index
h = L*sin(pi*n/L)./(pi*n); % ideal lowpass filter
h(33) = L;                % fill in the zero divison sample
h = hamming(65)'.*h;      % apply weighting window on h[n]

% interpolate:
y = filter(h,1,xe);     % y[n] is the interpolated signal