sampling - How to calculate FFT of a non-uniformly sampled signal?

Is it possible to calculate the FFT of a signal that was not sampled uniformly? I have an MPU6050 connected to an Arduino uno. I write the samples at the default frequency, approximately 20 milliseconds, but it's not always constant (see the image). Also see this related question.

This is my first time I approach to the FFT, and I have a civil engineer background, so is my first time with signal analysis, never studied Signal Theory before.

My interest for the FFT is to define the best low pass filter (example apply a butterworth but I don’t know how to choose the filter order and cutoff frequency.

Thanks for your time, thanks for you patience. I tried to solve this with the code below, but I think is only rubbish.

M=csvread('DATALOG.TXT'); 
n=length(M);

cleaning=500;
n=(n-cleaning);

MatFilt= M(cleaning:n,:);
m=length(MatFilt);
A= MatFilt(2:m,1);
B=MatFilt(1:m-1,1);
C=A-B;
meanrate= mean(C);
f= 1000/meanrate;

T= sum(MatFilt);    
TemporalSample= MatFilt(:,1);
Ax=MatFilt(:,2);
Axs=Ax*19.6/32768; %I converted the data in m/s^2
SignalLength= (MatFilt(m,1)-MatFilt(1,1))/1000;
plot(TemporalSample,Axs);
pause
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)');
NFFT = 2^nextpow2(m); % Next power of 2 from length of y

Y = fft(Axs,NFFT)/m;
f = f/2*linspace(0,1,NFFT/2+1);

Answer

You are essentially referring to nonuniform signal sampling in which data samples are not acquired at exact periods but have same random jitter deviations from their exact uniform timing positions.

Now under some suitable conditions (such as the jitter length being less than a period), there are very effective algorithms which can perfectly convert the nonuniformly sampled data into their equivalent uniform samples. (i.e., given the values of nonuniform samples they give a new set of values which correspond to the perfect uniform samples).

After this nonuiform to uniform conversion, you can then apply the usual FFT, which assumes that the data to be transformed was uniformly sampled.

However, before such a DSP attempt, I would suggest you understand the sampling and read out mechanism of the Arduino Uno and MPU 6050 6DOF sensor combination. That sensor can be read in either a single data read mode or in a burst buffer data read mode. In this second mode, I think it will provide a set of buffered data which are sampled by the MPU6050 unit itself as uniform as possible, compared to a sigle data read mode mastered by the Arduino.

From your data, there seems to be a drift of a few milliseconds per sample. I assume a 20 milli second of sampling period is quite achievable with less than 1 ms of timing accuracy using the Arduino in a proper configuration. So may be your main problem is on the software configuration?

EDIT: Try the following MATLAB code as is: NOTE: The following code MAY result in WORSE results than available before processing... The success highly depends on the application.

% NONUNIFORM SIGNAL DEFINITION
% ----------------------------
N=100;       % Converts a block of length N nonuniform samples to uniform
K=N;        % 
M = N*K;    % 
Tsam = 0.020; % Sampling period in seconds,

% tp represents the nonuniform sampling times
tp = [0:Tsam:(K-1)*Tsam] + Tsam*(rand(1,K)-0.5);

xnunds = sin(2*pi*1*tp);        % A simulated nonuniformly sampled data...

% GENERATE THE MATRIX AND ITS INVERSE
% -----------------------------------
A=zeros(K,K);       % initial matrix to be inverted
for p=1:K
    for q=1:K

        A(p,q)=sinc((tp(p)-tp(q))/Tsam);


    end
end
AI=inv(A);          % Inverse matrix is here


% APPLY THE CONVERSION FORMULA
% -----------------------------
xint=zeros(1,K);    % interpolate to uniform samples

for n=1:N:M           


    ACS=0;

    for m=1:K       % per nonuniform sample loop

        PSI=0;

        for q=1:K   % calculation of PSI for each m

            PSI = PSI + AI(q,m) * sinc( ((n-1)*Tsam/N - tp(q))/Tsam);


        end%q

        ACS=ACS+xnunds(m)*PSI;  % acumulation of multiplicant

    end%m

    xint( (n-1)/N + 1)=ACS;    % result for each n is assigned
end%n


figure,stem(xnunds);
hold on
stem(sin(2*pi*1*[0:K-1]*Tsam),'g');
title('nonuniform samples imposed on true samples');

figure,stem(xint,'r');
hold on;
stem(sin(2*pi*1*[0:K-1]*Tsam),'g');
title('interpolated uniform samples imposed on true uniform samples)