Register to reply 
Fourier Transform for unevenly sampled date 
Share this thread: 
#1
Feb1413, 10:53 AM

P: 5

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 xaxis (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 


#2
Feb1413, 12:00 PM

P: 534




#3
Feb1813, 06:57 AM

P: 5

Thanks for the reply. Perhaps I was not clear about a few aspects.
I need Y as a continous function of X where the the function consists of only sines andor 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 easytounderstand and easytoimplement nonuniform 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). 


#4
Feb1813, 08:04 AM

Engineering
Sci Advisor
HW Helper
Thanks
P: 7,293

Fourier Transform for unevenly sampled date
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 rethink whether this is a good way to approximate the data. 


#5
Feb1813, 10:01 AM

P: 5

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.



#6
Feb1813, 10:27 AM

Sci Advisor
PF Gold
P: 2,262

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 nonuniform 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. 


#7
Feb1813, 10:48 AM

Engineering
Sci Advisor
HW Helper
Thanks
P: 7,293

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


#8
Feb1913, 03:16 AM

P: 5




#9
Feb1913, 03:18 AM

P: 5

BTW, I found this (nonequispaced fast Fourier transform or nonuniform fast Fourier transform): http://wwwuser.tuchemnitz.de/~potts/nfft/



#10
Feb1913, 08:39 PM

P: 699

If it were me I would do a leastsquares 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 leastsquares 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! % leastsquares 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 (leastsquares). 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 


#11
Feb1913, 09:26 PM

P: 589

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. 


Register to reply 
Related Discussions  
Computing the Hilbert transform via Fourier transform  Math & Science Software  2  
Fourier transform, Fourer Integral transform  Calculus & Beyond Homework  5  
Discrete Fourier transform of sampled continuous signal  Engineering, Comp Sci, & Technology Homework  0  
Fourier Series / Fourier Transform Question  Electrical Engineering  6  
The difference between Fourier Series, Fourier Transform and Laplace Transform  General Physics  1 