time integration via frequency domain

80 views (last 30 days)
Hekseli
Hekseli on 25 Jun 2014
Commented: Antonio Leanza on 13 Mar 2023
Hey!
I have been trying to get this integration by fft and dividing by 2*pi*f*i working but the amplitude is wrong all the time. Basically my code is the same as in this example: http://www.mathworks.com/matlabcentral/answers/17611-how-to-convert-a-accelerometer-data-to-displacements
So
function int_time_data = integrate(data,dt)
N1 = length(data);
N = 2^nextpow2(N1);
if N > N1
data(N1+1:N) = 0; % pad array with 0's
end
df = 1 / (N*dt); % frequency increment
Nyq = 1 / (2*dt); % Nyquist frequency
freq_data = (fftshift(fft(data)));
int_freq_data = zeros(size(freq_data));
f = -Nyq:df:Nyq-df;
for i = 1 : N
if f(i) ~= 0
int_freq_data(i) = freq_data(i)/(2*pi*f(i)*sqrt(-1));
else
int_freq_data(i) = 0;
end
end
int_time_data = ifft(ifftshift((int_freq_data)));
int_time_data = int_time_data(1:N1);
end
dt = 0.01; % seconds per sample
N = 512; % number of samples
t = 0 : dt : (N-1)*dt; % in seconds
wave_freq = 17.1; % in Hertz
data = sin(2*pi*wave_freq*t);
integrated_time_data = integrate(data,dt);
integrated_time_data_analytical = -1/(2*pi*wave_freq)*cos(2*pi*wave_freq*t);
plot(t,integrated_time_data);
hold on;
plot(t,integrated_time_data_analytical,'-r');
But this produces wrong results. What is wrong?
  1 Comment
Antonio Leanza
Antonio Leanza on 13 Mar 2023
In this code the detrendization via high pass filtering is missing.

Sign in to comment.

Answers (3)

Star Strider
Star Strider on 25 Jun 2014
You did not say what was wrong about the amplitude. I did not run that code, but see if changing the freq_data line to:
freq_data = (fftshift(fft(data)/length(data)));
solves the problem.
  2 Comments
Star Strider
Star Strider on 25 Jun 2014
That may be required, since according to the ifft documentation, it also normalizes by dividing the computed ifft by the length of the argument.
Antonio Leanza
Antonio Leanza on 13 Mar 2023
Edited: Antonio Leanza on 13 Mar 2023
The original code works well regarding that line.

Sign in to comment.


Hekseli
Hekseli on 25 Jun 2014
Thanks for the hint. It didn't solve the problem totally. When dividing by the length in fft I assume I should multiply by the length after the ifft? Otherwise the amplitude is far too small.
Now after the multiplication I get the following result which is attached

Michael
Michael on 27 Aug 2014
Hekseli,
You sometimes need to perform a high pass filter on your data when you convert it from acceleration to position. That should get rid of any low frequency drift you see in your results.

Categories

Find more on Fourier Analysis and Filtering in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!