Zero-padding in the time domain corresponds to interpolation in the Fourier domain. It is frequently used in audio, for example for picking peaks in sinusoidal analysis.
While it doesn't increase the resolution, which really has to do with the window shape and length. As mentioned by @svenkatr, taking the transform of a signal that's not periodic in the DFT size is like multiplying with a rectangular window, equivalent in turn to convolving it's spectrum with the transform of the rectangle function (a sinc), which has high energy in sidelobes (off-center frequencies), making the true sinusoidal peaks harder to find. This is known as spectral leakage.
But I disagree with @svenkatr that zero-padding is causing the rectangular windowing, they are separate issues. If you multiply your non-periodic signal by a suitable window (like the Hann or Hamming) that has appropriate length to have the frequency resolution that you need and then zero-pad for interpolation in frequency, things should work out just fine.
By the way, zero-padding is not the only interpolation method that can be used. For example, in estimating parameters of sinusoidal peaks (amplitude, phase, frequency) in the DFT, local quadratic interpolation (take 3 points around a peak and fit a parabola) can be used because it is more computationally efficient than padding to the exact frequency resolution that you want (would mean a much larger DFT size).
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 SebastianR@gmx.de if you really want to implement this.
Best Answer
Sure, you can use a radix-2 FFT to compute FFTs for lengths not a power of 2 (but it is not as efficient as using methods specifically tailored to the factors of the sequence length). However, it's not as easy as merely padding the original array. The right way of going about it, due to Rabiner, Schafer, and Rader, is termed the "chirp-z transform". (This is a slightly condensed version of the paper previously linked to.)
To motivate the chirp-z idea, consider the DFT
$$X_k=\sum_{n=0}^{N-1} x_n \left(e^{-2\pi i/N}\right)^{k n},\qquad k = 0,\dots,N-1$$
The key step, attributed to Bluestein, considers the algebraic identity
$$kn=\frac12(k^2+n^2-(k-n)^2)$$
Substituting into the DFT and expanding out the powers yields
$$\begin{align*}&\sum_{n=0}^{N-1} x_n \left(e^{-2\pi i/N}\right)^{\frac12(k^2+n^2-(k-n)^2)}\\=&\left(\left(e^{-2\pi i/N}\right)^{k^2/2}\right)\sum_{n=0}^{N-1} \left(x_n \left(e^{-2\pi i/N}\right)^{n^2/2}\left(e^{-2\pi i/N}\right)^{-\frac{(k-n)^2}2}\right)\end{align*}$$
and we then recognize that the summation expression is precisely in the form of a convolution. The factor $\left(e^{-2\pi i/N}\right)^{k^2/2}$ is what is termed as a "chirp"; hence the name "chirp-z transform".
Chirp-z thus consists of three stages:
MATLAB's Signal Processing toolbox implements this algorithm as the function
czt(). See the code in the corresponding M-file for implementation details.The Code
I (the original poster) ended up making a Python version using SciPy, so I thought I'd edit this answer to include the code: