I have a waveform which I would like to take the fourier transform of. However, the simulator which generated the waveform uses an adaptive algorithm to determine the timestep for each calculation.
Now, I figure I can regularize the timestep to first order by creating a new timebase and for each point in the new base calculating a linear interpolation from the two nearest points in the original sequence, but this seems inelegant, and potentially incorrect.
I also considered expanding the original sequence with whatever the smallest time-step is, replicating data points to fill-in the missing data (effectively creating steps), and then filtering the resulting fourier transform above some frequency (or simply ignoring data above some frequency).
Either of these will probably work with the data I have, since it is heavily oversampled for the frequency range I am interested in, but I worry that both would cause distortions if that were not the case.
Is there a better method?
Best Answer
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.