There may be literature on this, and I'm no expert, but I believe there is a simple method to directly compute what you want:
First, FFT is just a fast algorithm for the DFT. For the following method, you will need to compute the DFT coefficients directly. Recall the formula for a single coefficient (slightly modified from Wikipedia to make it match what I say below):
$$X_k := \frac{1}{N} \sum_{n=0}^{N-1} x_n e^{-\pi i k \frac{2n+1}{N}}$$
You can interpret this as a discretization of the integral (hope I'm getting this right):
$$Y_k := \int_0^1 y(t) e^{-2 \pi i k t} dt$$
While $x$ is an $N$-tuple, $y$ is supposed to be the continuous function on $[0,1]$ that $x$ approximates. Each $n=0,\dots,N-1$ is associated with the interval $[a_n,b_n]:=[\frac{n}{N},\frac{n+1}{N}]$, so the union of all intervals is $[0,1]$. It also makes sense to assume there is a specific $t_n\in[0,1]$ such that $x_n = y(t_n)$, say $t_n := \frac{a_n + b_n}{2}$. (This is what I implicitly used in the first formula, making it a midpoint sum approximation. Please ask if this is not clear.)
Now, the input you have is basically a different discretization $z$ of $y$, which is an $M$-tuple ($M$ possibly being different from $N$), and where $a_m$, $b_m$, and $t_m$ for $m=0,\dots,M-1$ are given by your variable timestep. (You have to decide for yourself whether the interval boundaries or the measurement points are given; however, you have to make sure that $a_0=0$ and $b_{M-1}=1$.) What you need is a sum that appropriately discretizes the above integral, which would (hopefully, again) be:
$$Z_k := \sum_{m=0}^{M-1} (b_m-a_m) z_m e^{-2 \pi i k t_m}$$
If the timestep is completely variable, you can get nonzero coefficients for arbitrarily large $k$, but you can stop calculating wherever you want, of course.
Now, computing these sums is not nearly as fast as an FFT, but there are several things you can do to make it computationally feasible. The most important thing is to calculate the sines and cosines just once for each $m$, i.e. iterate over $m$ at the outer level, not $k$. You should pre-calculate sine/cosine tables if you have a fixed "time resolution" that every timestep is a multiple of; otherwise use sine/cosine implementations that favor speed over accuracy.
If you want overlap between consecutive DFT "windows" (not sure if this is the correct word here), there is a way to reuse the part of each sum that is in both windows. Another way to reduce the time overhead is to use different window lengths for each coefficient, which gives you a variable (e.g. logarithmical) frequency resolution. I don't have time to explain this now, but you can mail me for details at [email protected] if you really want to implement this.
