I know that if we have some data representing some wave, for example image line values, we can use fourier transform to get frequency function of that wave. But we have N values at points x=0...N1 And we get only N frequencies at the output. So I want to analyze the wave everywhere in the range [0, N1] For example at point u = 1.5. How can I do it?
You need to have some a priori knowledge about your wave between and beyond the sampling points to get a meaningful guess about the full Fourier transform. The $N$ values that you get doing the discrete Fourier transform have anything to do with the continuous Fourier transform only for indices much less than $N$. Note that different assumptions will lead to different answers for large frequences. If you need only relatively low frequences, your signal is compactly supported in time, and your sampling points cover the entire support, you can safely interpolate the values of the discrete FT to guess the continuous one but there is no way to get reliable high precision values for the continuous Fourier transform at frequences comparable to $N$. The $N$point resolution on $[0,1]$ is just not high enough to catch those frequences without an error. Also note that you want the frequences not on "on the interval $[0,N1]$" as you wrote, but rather on the interval $[M,M]$ where $M\ll N$. For decent signals, the frequences with indices close to $N$ in the discrete $FT$ actually catch the negative low frequences in the continuous FT, not the high ones. 


Strictly speaking, you can't. That is, you have to interpolate to do it, and to do that in a meaningful way requires some knowledge of the problem domain. However, the easiest way to do so is to note that the given function is naturally written in terms of its discrete Fourier transform as a trigonometric polynomial, so you can just go ahead and evaluate that trigonometric polynomial at any intermediate point. The same thing goes with the function and its Fourier transform interchanged, of course. Though it does not really make good mathematical sense to do this interpolation in both domains simultaneously. 


Here is one answer from the electrical engineering(or, image processing) point of view. However this answer much nonmathematical and more for practical illustration. Since this is a nonmathematical answer, I have made it community wiki so that I don't get reputation. In digital signal processing, you sample an analog signal into digital encoding, and let your microprocessor work with it. When you do sampling, you clearly lose some information and there is no way you can recover it. When you become too greedy and actually try to get more, the problem of "aliasing" will happen. Here is the wikipedia link for aliasing. For image processing, situation is even better for an illustration, since you can actually see wavy patterns. These are called Moiré patterns, and you must have seen it sometime in the real world. Here is a wikipedia link. 


If what you need is a simple practical method to do the interpolation,
then just multiply the "time domain" samples by the linearphase
signal $\{\exp\left(\alpha\cdot i 2\pi k /N\right)\}_{k=0}^{N1}$, where
$\alpha\in(0,1)$ is the subsample 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
continuoustime signal $\{X(t)\}_{t\in \mathbb{R}}$, as in
$x_k=X(kT_s)$ for $k=0,\ldots, N1$, where $T_s\in \mathbb{R}_{++}$ is
the sampling period. If the continuoustime signal has a compactly 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 

