Resampling scattered data with MATLAB's $\tt interp$ and $\tt resample$

Question

I am recording acceleration data with an MPU6050 connected to a Arduino1 and stored on a SD. Here you can find the code.

I need to calculate the FFT of an acceleration signal that was not sampled uniformly, so i have to resample my signal. After this suggestion I tried to study how to resample the signal.

From Mathworks documentation here it seems that I have to know the frequency of my signal and a nominal frequency. I read and studied the function resample, but I don’t understand how it works very well. In my specific case how to define

x= sin(2*pi*f*T) 

But in my code below I tried to implement the function (to give you more information for a detailed answer and prove my efforts, but I don't have enough reputation to post all the graph).

Then I studied gettdata “Scattered interpolant class” finding not very useful in my situation, they seems to work for 2 or more dimensions.

So I followed another way, I used the function interp1. I defined my equispaced vector using the mode and median value of my sample frequency, I discarded the average to avoid the influence of any peaks.

• Is a correct approximation?
• What kind of errors can occur?

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.

I don't have enough reputation to post more than 2 links sorry for that.

filename= uigetfile ('.txt');
fileID = fopen (filename);
logmpu6050 =csvread(filename);
fclose (fileID); 
%Starting creating the specific Vectors
%Record Time in millisecond 
time=logmpu6050(:,1);
%The x y z are converted to m/s^2
confactor=19.6/32768;
ax=logmpu6050(:,2);
confactor=19.6/32768;
ax=ax*confactor;
ay=logmpu6050(:,3);
ay=ay*confactor;
az=logmpu6050(:,4);
az=az*confactor;
%Define the sample rate subtracting from Sampletime i+1 Sampletime i
n=length(time);
for i = 1:n-1
samplerate(i) = time(i+1)-time(i);
end;

%I try to define the best frequency for "resample" function based on mode and median value 
%I will not use average to avoid conditioning due to extreme value 
f1=mode(samplerate);
f2=median(samplerate);
x1= sin(2*pi*f1*time);
x2= sin(2*pi*f2*time);
plot(time,x1,'.');
figure
plot(time,x2,'.');  

%resampling at time mode record frequency
[test1,test2]= resample(az,time,f1);
%resampling at time median record frequency
[test3,test4]= resample(az,time,f2);
figure
plot(test2, test1)
figure
plot(test4,test3)     

f1=int32(mode(samplerate));
f2=int32(median(samplerate));

%Try to implement the matlab function "interp1"
%My input are 
%x-> time when the value is recorded "time"
%v-> the acceleration value "az" 
%xq-> I tried two define two type of equispaced vectors
%xq1-> Defined by the total recording time divided by the mode
%xq2-> Defined by the total recording time divided by the median 
totrectime= int32(logmpu6050(n,1)-logmpu6050(1,1)) %milliseconds
numsamples1= idivide(totrectime,f1); %Rounded toward zero to the nearest integers.
numsamples2= idivide(totrectime,f2); 
f1double=double(f1);
f2double=double(f2);
for i = 1: numsamples1
xq1(i)= f1double*i;
end;
for i = 1:numsamples2
xq2(i)= f2double*i;
end;
xq1double=double(xq1);
xq2double=double(xq2);

vq1 = interp1(time,az,xq1double);
vq2 = interp1(time,az,xq2double);
figure
plot(xq1double,vq1);
figure
plot(xq2double,vq2);

Answer 1

In your question, you state that your ultimate goal is to compute the discrete Fourier transform (DFT) from nonuniformly-sampled input data. It is not necessary to resample your data at all. In fact, there are several fast algorithms, known variously as the unequally spaced fast Fourier transform (USFFT) and the nonuniform fast Fourier transform (NUFFT) that do exactly what you want: take nonuniform samples as input and produce equally spaced samples of the DFT as output. The remainder of your question, dealing with issues related to MATLAB's resampling routines, are a red herring.

Here are some publications describing the USFFT algorithms:

1. G. Beylkin, On the Fast Fourier Transform of Functions With Singularities, Applied and Computational Harmonic Analysis, 2, pp. 363-381, 1995. Link to preprint.

2. A. Dutt and V. Rokhlin, Fast Fourier Transforms for Nonequispaced Data, SIAM J. Sci. Comput., 14, pp. 1368, 1993.

There are many commercial and free implementations of these algorithms. You can access some of them here.

EDIT: Here is a "quick and dirty" approach that may help you.

Based on your problem description and comments, you are just trying to get a sense of the spectral content of your signal, which does not require high accuracy. Here is one way to formulate your problem so that the spectrum can be easily estimated using any of the USFFT implementations in the above link.

You have a signal $f(t)$ defined on the $T$-duration interval $t\in[-T/2,T/2]$, and you wish to estimate its Fourier transform $$\hat{f}(\omega)=\int_{-T/2}^{T/2}f(t)e^{-i 2 \pi \omega t}\,dt$$ at a number of values of $\omega$. You have collected $N$ samples of your function $$f_n=f(t_n),\quad n=1,\ldots,N$$ at the unequally spaced time locations $t_n$, $n=1,\ldots,N$.

In general, you would solve this by formulating a least squares problem for the Fourier transform $\hat{f}$, and using the USFFT to rapidly solve the system of equations, but since you do not require high accuracy, here is a "quick and dirty" solution.

Use the trapezoid rule from Calculus 101 to discretize the above integral, giving the approximation $$\hat{f}(\omega) \approx \frac{\Delta_1}{2}f_1 e^{-i2\pi\omega t_1}+\left(\sum_{n=2}^{N-1}\frac{\Delta_{n-1}+\Delta_{n}}{2}f_n e^{-i2\pi\omega t_n}\right)+\frac{\Delta_{N-1}}{2}f_N e^{-i2\pi\omega t_N}$$ where $$\Delta_{n}=t_{n+1}-t_n$$ is the temporal spacing between samples $n+1$ and $n$.

Define the new sequence $$g_n=\begin{cases}\frac{\Delta_1}{2}f_1, & n=1 \\ \frac{\Delta_{n-1}+\Delta_{n}}{2}f_n, & n=2,\ldots,N-1 \\ \frac{\Delta_{N-1}}{2}f_N, & n=N\end{cases}$$ Now pick some frequency grid you care about. For example, say you want $2000$ equally spaced frequencies between $0$ Hz and $25$ Hz, and define $$\omega_m = \frac{m-1}{2000-1}25,\quad m=1,\ldots,2000.$$

Now your problem amounts to evaluating the sum $$\hat{f}(\omega_m) = \sum_{n=1}^N g_n e^{-i2\pi\omega_m t_n},\quad m=1,\ldots,M.$$

This looks like an FFT, except the $t$-values are unequally spaced and the $\omega$-values are unusual. If you do not care about speed, then just evaluate the sum directly using your favorite programming language. It will take a minute or two of computer time, but so what? If you are just doing this a few times to get a sense of what your signal's spectrum looks like, then it does not matter how fast it is. If at some point in the future you decide that you do care about speed, then the sum can be evaluated extremely rapidly using a single call to a USFFT library, which you can download from the link I posted (or you can use one of the many commercial implementations). Making the switch to the USFFT will only require you to modify a handful of lines of code.