Fourier Transform for unevenly sampled date

[bendanish]

Dear people,

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.

In other words, in the x-axis (angle of rotation in degrees for three rotations) I have data which looks like this:
x = [0 55 125 190 240 300 350 65 120 185 240 355 50 110 180 235 355 0];
y = [0 1.2 2 2.5 2.1 1.1 0.1 0.9 1.9 2.6 2 1.0 0 0.8 1.9 2.5 1 0];

I basically need the Fourier transform of Y.

Also, its not possible for me to do any sort of interpolation to make X look "nice".

Would really appreciate if anyone can help me out with this, either my recommending some books or papers that deal with these sort of problems or by providing a code that can do this.

Thanks.

Danish

 

[milesyoung]

I measure the angle of rotation and the verticle displacement at a steady rate of 10,000 samples per second.

You're simultaneously sampling angle and displacement at 10 kHz? If so, then you already know how the vertical displacement varies with the angle of rotation.

I basically need the Fourier transform of Y.

Why is it you think you need to apply the Fourier transform?

 

[bendanish]

Thanks for the reply. Perhaps I was not clear about a few aspects.

I need Y as a continuous 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).

 

[AlephZero]

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.

 

[bendanish]

Yes, I don't need all 6000 coefficients. Unfortunately, I have to get the Fourier coefficients because my boss told me to:( Otherwise I would also use other interpolation techniques.

 

[f95toli]

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.

 

[AlephZero]

Unfortunately, I have to get the Fourier coefficients because my boss told me to:( Otherwise I would also use other interpolation techniques.

If your boss told you to jump out of a window on the 20th floor, would you do that as well?

If all else fails, just bounce the question back and ask him/her how to do the calcs

 

[bendanish]

If your boss told you to jump out of a window on the 20th floor, would you do that as well?

Yes, its exactly the same thing.

 

[bendanish]

BTW, I found this (nonequispaced fast Fourier transform or nonuniform fast Fourier transform): http://www-user.tu-chemnitz.de/~potts/nfft/

 

[jasonRF]

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!

% least-squares fit
%
% model:
% let y(t) = a0 + a1*cos(t) + a2*cos(2*t) + ... + aM*cos(M*t) + b1*sin(t) + b2*sin(2*t) + ... + bM*sin(M*t)
%
% have values of t (angle) and y.

% data from post
t = [0, 55, 125, 190, 240, 300, 350, 65, 120, 185, 240, 355, 50, 110, 180, 235, 355, 0].'; %angle
y = [0, 1.2, 2, 2.5, 2.1, 1.1, 0.1, 0.9, 1.9, 2.6, 2, 1.0, 0, 0.8, 1.9, 2.5, 1, 0].'; %displacement

N = length(t); %number of points

% select maximum "frequency"
M = 4;

if (2*M+1>N),
disp('not enough data for this fit!')
end;

% 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);

% first column of A corresponds to a0 = const
A(:,1) = ones(N,1);

% build next M colums corresponding to cosine terms
for ii=1:M,
A(:,ii+1) = cos(ii*t*pi/180);
end

% build next M colums corresponding to sine terms
for ii=1:M,
A(:,M+ii+1) = sin(ii*t*pi/180);
end% now do least squares - use Octave (or matlab) backslash for best numerics
c = A\y;

% break out coeffs
a0 = c(1); %constant term
a = c(2:M+1); %a(1) = a1, a(2) = a2, etc.
b = c(M+2:2*M+1); %b(1) = b1, b(2) = b2, etc.% now can use this model to evaluate the vertical displacement y for any theta. For example
theta = [0:1:359].'; %angles at which to evaluate model
Nt = length(theta);
ymodel = zeros(Nt,1);

ymodel(:) = a0; %constant term
for ii=1:M,
ymodel = ymodel + a(ii)*cos(ii*theta*pi/180);
ymodel = ymodel + b(ii)*sin(ii*theta*pi/180);
end

% plot model as well as data
figure
plot(theta, ymodel); %model using line
hold on;
plot(t, y, '*'); %data using markers

 

[the_emi_guy]

You may want to consider unwrapping the angle. In other words, don't have it roll back to 0 after 360.

For example, have it go up to 4320 after 12 rotations and take FFT of first 4096 samples.