time integration via frequency domain

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?

[Antonio Leanza]

In this code the detrendization via high pass filtering is missing.

Answer 1

[Star Strider]

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.

[Star Strider]

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]

The original code works well regarding that line.

Answer 2

[Hekseli]

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

Answer 3

[Michael]

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.