Fourier Transform for unevenly sampled date

I am trying to analyze data from test bench which consists of a magnetically levitated spindle. We have a rotor/spindle which rotates and moves vertically up and down as it rotates. I measure the angle of rotation and the verticle displacement at a steady rate of 10,000 samples per second.

The issue is that the rotational speed is not perfectly constant so for one rotation of the of the rotor I somtimes get 5k samples, sometimes 6k samples.

I need to express the vertical movement as a function of the rotation angle.

I need Y as a continous function of X where the the function consists of only sines and-or cosines. So basically I need the fourier coefficients.

I need this in order to be able to approximate high derivatives of the the function Y. I will later compare the results with some numerical methods of derivate approximation.

In essence, I am looking for an easy-to-understand and easy-to-implement non-uniform DFT or FFT. I record 6000 data points per revolution so DFT might be too cumbersome (N*N computations means 36 million in this case).

Presumably you are not trying to find all 5000 or 6000 fourier coefficients, but only the lowest frequency ones (i.e. at most a few hundred).

You can do this as a least squares approximation problem to fit
##y(x) = a_1 f_1(x) + a_2 f_2(x) + \dots## where the ##f_i## just happen to be sine and cosine functions.

For equally spaced points this is exactly equivalent to the standard DFT (but of course the FFT algorithm is faster).

Solve the least squares problem with an efficient method (e.g. based on the QR algorithm or singular value decomposition) rather than just forming the normal equations and solving them.

For unequally spaced points, note that the coefficients will depend on the number of terms you include in the fit (i.e. if you do two calculations to find 50 coefficents and 100 coefficients, the first 50 will not be exactly the same both times). But if the differences are significant, you probably need to re-think whether this is a good way to approximate the data.

AFAIK the only way to do this if you want to use an FFT is to first interpolate the data and then apply a standard FFT, evenly spaced data is an underlying assumption for all FFT algorithms and the usual way to achieve this with non-uniform data is to use simple interpolation.
Hence, if you also exclude a full DFT there is no "simple" way of doing it and you have to start considering more clever algorithms for calculating the periodogram, but the problem with that is that the results are not always easy to interpret.
That is, you need to find a good DSP book.

If it were me I would do a least-squares fit to get a nice model. Presumably the vertical displacement should be a function of only the angle, as opposed to also depending upon the history of the dynamics?

Anyway, assuming that y is only a function of the angle, you can do a least-squares fit. Use sinusoids since the model must be periodic with period of 360 degrees. Note that this approach may force you to make a HUGE matrix. But with modern computers you will likely be okay.

EDIT: note that if you had evenly sampled data (and only one value for a given angle) then the sinusoids are orthogonal, which means that you do not have to form the giant A matrix below since each coefficient is independent! But with unevenly sampled data you have no such luck!

It was easier to write Octave (matlab) code than write latex. Hopefully it is self explanatory. Ask questions!

% build matrix. Do the slow way using loops to show how this works.
%
% want to solve A*c = y, where y is the measured data, c are the coeffs of fit (a0,a1, ..., aM, b1, ..., bM)
% since in general is an overdetermined and inconsistent system, we will find the solution
% that minimizes the square error (least-squares).
A = zeros(N, 2*M+1);