If what you need is a simple practical method to do the interpolation, then just multiply the "time domain" samples by the linear-phase signal $\{\exp\left(\alpha\cdot i 2\pi k /N\right)\}_{k=0}^{N-1}$, where $\alpha\in(0,1)$ is the sub-sample shift in the "frequency domain." (I'm not sure if an additional constant of absolute value 1 is required here). As already noted in previous answers, the motivation for this is the assumption that your samples vector $\{x_k\}\\\$ comes from sampling some continuous-time signal ${X(t)}{t\in \mathbb{R}}$, as in $x_k=X(kT_s)$ for $k=0,\ldots, N-1$, where $T_s\in \mathbb{R}{++}$ is the sampling period. If the continuous-time signal has a compactly supported Fourier transform, than by Shannon's sampling theorem (http://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem) $X$ is determined by the infinite sequence ${X(kT_s)}_{k\in \mathbb{Z}}$ for small enough $T_s$. Since we only have $N$ entries of the infinite sequence, we loose some information on $X$. But if $X$ decays rapidly (rather than having a compact support), then at least intuitively we don't loose much (I know that this is a dangerous and imprecise statement :) ).
For example, try the following Matlab lines:
ttt=-32:31;
x_time=exp(-ttt.^2/10);
figure;plot(x_time);
x_freq = fft(x_time);
figure;plot(fftshift(abs(x_freq)));
figure
for alpha=-4:0.5:4
x_time_2 = x_time .* exp(2 * pi * i * (0:63)/64 * alpha);
x_freq=fft(x_time_2);
plot(fftshift(abs(x_freq)));
axis([26,37,0,10]);
grid on;
pause(1);
end