I am trying to obtain velocity and position data using the Fast Fourier Transform. The algorithm I use is described in the below figure; enter image description here
I have achieved to succesfully calculate it offline (when the accelaration measurement data is gathered over a time interval and the calculation of velocity and position for the same interval is done after all the data is gathered.). The function I use is as below;
function [velocity,position]=acc2velpos(acc,dt,f_low)
nsamples = length(acc(:,1));
domega = 2*pi/(dt*nsamples);
acchat = fft(acc);
% Make frequency array:
Omega = zeros(nsamples,1);
% Create Omega:
if rem(nsamples,2) == 0
nh = nsamples/2;
Omega(1:nh+1) = domega*(0:nh);
Omega(nh+2:nsamples) = domega*(-nh+1:-1);
else
nh = (nsamples-1)/2;
Omega(1:nh+1) = domega*(0:nh) ;
Omega(nh+2:nsamples) = domega*(-nh:-1);
end
% High-pass filter:
n_low=floor(2*pi*f_low/domega);
acchat(1:n_low)=0;
acchat(nsamples-n_low+1:nsamples)=0;
% Multiply by omega^2:
% Skip the measurement in the first time step since it might be inaccurate
poshat(2:nsamples) = -acchat(2:nsamples) ./ Omega(2:nsamples).^2;
velhat(2:nsamples) = acchat(2:nsamples) ./ (Omega(2:nsamples)*1j);
% Inverse Fourier transform:
velocity = real(ifft(velhat));
position=real(ifft(poshat));
% --- End of function ---
end
I used the acceleration input as below, sampling frequency Fs = 256 Hz, sample size N = 2560 (sample length 10 s), and cut-off frequency 1 Hz
Using this dataset I obtained the same results with the article (Gu Mingkun,Lü Zhenhua,Identification of a Mechanism's Vibration Velocity and Displacement Based on the Acceleration Measurement 2010) proposed this algorithm over the 10 second time interval. enter image description here
However my main objective is to use this algorithm to calculate velocity and displacement online. What I did to adapt it is to use a Simulink model as below; enter image description here
The Matlab Function block code is below;
function [vel,pos, veldumcheck]=acc2velpossimulink(acc)
dt=10/2560;
f_low=1;
%velstore=zeros(2560,2560);
persistent inputArray;
% Initialize the array on the first function call
if isempty(inputArray)
inputArray = zeros(1,1);
end
% Store the input in the array
inputArray = [inputArray; acc];
if length(inputArray)<1000
vel=0;
pos=0;
veldumcheck=0;
else
nsamples = length(inputArray)-1;
domega = 2*pi/(dt*2560);
acchat = fft(inputArray(2:end));
% Make frequency array:
Omega = zeros(nsamples,1);
% Create Omega:
if rem(nsamples,2) == 0
nh = nsamples/2;
Omega(1:nh+1) = domega*(0:nh);
Omega(nh+2:nsamples) = domega*(-nh+1:-1);
else
nh = (nsamples-1)/2;
Omega(1:nh+1) = domega*(0:nh) ;
Omega(nh+2:nsamples) = domega*(-nh:-1);
end
% High-pass filter:
n_low=floor(2*pi*f_low/domega);
acchat(1:n_low)=0;
acchat(nsamples-n_low+1:nsamples)=0;
velhat=zeros(1,nsamples);
velhat = complex(velhat, zeros(size(velhat), class(velhat)) );
poshat =zeros(1,nsamples);
poshat = complex(poshat, zeros(size(poshat), class(poshat)) );
% Multiply by omega^2:
% Skip the measurement in the first time step since it might be inaccurate
poshat(2:nsamples) = -acchat(2:nsamples) ./ Omega(2:nsamples).^2;
velhat(2:nsamples) = acchat(2:nsamples) ./ (Omega(2:nsamples)*1j);
% Inverse Fourier transform for velocity:
veldum=zeros(1,nsamples);
veldum(:) = real(ifft((velhat)));
veldumcheck=veldum(:);
vel=veldum(end);
% Inverse Fourier transform for position:
posdum=zeros(1,nsamples);
posdum(:) = real(ifft(poshat));
pos=posdum(end);
% --- End of function ---
end
end
The main difference is it takes the acceleration data at each time step, stores in an array and calculates the velocity until that time interval and lastly outputs the velocity data at last time step. However the result is very different than what I calculated offline.
I figured out the FFT algorithm I described above tend to result with a negative value at the last time step. For example when I run the Matlab code for 10 second time interval I see that the velocity has a positive value between t=4.5s and t=4.85s however when I run the code for 4.75 second with the same input data set, velocity has a negative value at the last time steps of the time interval. And this is the same in most of the cases. As a result, in the Simulink function code outputing the last element of the velocity array result with a negative and inaccurate value.
Any idea on how to achieve online calculation of this problem is more than welcome.