How do different ways of interpolating complex numbers affect the discrete Fourier transform?

Question

I have an array of 1000 complex numbers $\vec b$. Here is an example of the norms of the numbers in $\vec b$:

Taking the fast Fourier transform of $\vec b$, then the norm, gives the frequency space spectrum:

I want to increase the range of this spectrum. Since I do not have more time domain data, I want to interpolate the values in $\vec b$. Obviously, the high frequency information will depend on the method of interpolation. Here is the previous frequency spectrum with double the amounts of points in $\vec b$, linearly interpolated:

Changing linear interpolation to spline, the frequency spectrum keeps about the same shape in the low frequencies (~800-1200) but changes slightly outside of that range, as expected. By assuming that the phase changes slowly between points, I used another interpolation scheme: $z_{interp} = e^{i\phi_{interp}}(|z_1|^2+|z_2|^2)/2$ where $\phi_{interp} = (arg(z_1)+arg(z_2))/2$, and this is the frequency spectrum that came out

I do not understand how this is possible as the low frequencies should not have been affected so drastically. Obviously there is something I do not understand about the relationship between interpolation and Fourier transforms.