4

I have time series with accelerometer data which I want to integrate to velocity and position time series. This is done using FFT, but the two methods in Matlab and Python give me different results.

Matlab Code:

nsamples = length(acc(:,1));
domega = 2*pi/(dt*nsamples);      
acchat = fft(acc);

% Make frequency array:
Omega = zeros(nsamples,1);

% Create Omega:
if even(nsamples)
   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:
acchat(2:nsamples) = -acchat(2:nsamples) ./ Omega(2:nsamples).^2;

% Inverse Fourier transform:
pos = real(ifft(acchat));

% --- End of function ---

Python Code:

import numpy as np


def acc2velpos(acc, dt):

    n = len(acc)

    freq = np.fft.fftfreq(n, d=dt)
    omegas = 2 * np.pi * freq
    omegas = omegas[1:]

    # Fast Fourier Transform of Acceleration
    accfft = np.array(np.fft.fft(acc, axis=0))

    # Integrating the Fast Fourier Transform
    velfft = -accfft[1:] / (omegas * 1j)
    posfft = accfft[1:] / ((omegas * 1j) ** 2)

    velfft = np.array([0] + list(velfft))
    posfft = np.array([0] + list(posfft))

    # Inverse Fast Fourier Transform back to time domain
    vel = np.real(np.fft.ifft(velfft, axis=0))
    pos = np.real(np.fft.ifft(posfft, axis=0))

    return vel, pos

The two codes should generally give the exact same results, but when I plot a comparison this is what I get

Acceleration to Velocity

enter image description here

Acceleration to Position

enter image description here

Any idea why the Python results are not identical to Matlab results in Position? It is crucial for me to have the same results as I am using these acceleration measurements from experiments to identify the motion of a barge.

I also have a second version in Python where I try to include the filter which is in the Matlab code. This also gives different results, much like the one without a filter in Python.

def acc2vel2(acc, dt, flow):

    nsamples = len(acc)
    domega = (2*np.pi) / (dt*nsamples)
    acchat = np.fft.fft(acc)

    # Make Freq. Array
    Omega = np.zeros(nsamples)

    # Create Omega:
    if nsamples % 2 == 0:
        nh = int(nsamples/2)
        Omega[:nh] = domega * np.array(range(0, nh))
        Omega[nh:] = domega * np.array(range(-nh-1, -1))

    else:
        nh = int((nsamples - 1)/2)
        Omega[:nh] = domega * np.array(range(0, nh))
        Omega[nh:] = domega * np.array(range(-nh-2, -1))

    # High-pass filter
    n_low = int(np.floor((2*np.pi*flow)/domega))
    acchat[:n_low - 1] = 0
    acchat[nsamples - n_low:] = 0

    # Divide by omega
    acchat[1:] = -acchat[1:] / Omega[1:]

    # Inverse FFT
    vel = np.imag(np.fft.ifft(acchat))

return vel

This still gives a little bit different than the Matlab code. Suggestions?

Improve this question

2 Answers 2

Reset to default
4

It looks like you have a high pass filter in the matlab code and not in the python code. Which makes sense given the difference between your python and matlab position results.

Your python position wave appears to oscillate up and down at a low frequency, indicating that some low frequency component in the frequency domain was not filtered out. The high pass filter in your matlab code removed the low frequency component.

Improve this answer
Sign up to request clarification or add additional context in comments.

4 Comments

I have another version of the Python Code where I have attempted to include the high-pass filter, but I still get a difference. This code is here.
Kakemonster, in your python code which includes the filter, are the results unchanged from when there was no filter or are they different from your original python results? Also, go through your code and look at the values for arrays and variables and see at what point the matlab and python values begin to differ.
Hey! I figured it out! It is like you said, the High Pass filter made the difference. After successfully implementing the same filter in python the results are exactly the same. Thanks for the help!!
No problem, glad to help. Would you please add your working code to your question as an edit? This will help other people with similar issues in the future.
0

Implementing the high-pass filter in the Python script solved the problem:

def acc2velpos(acc, dt, f_low):

    n = len(acc)

    freq = np.fft.fftfreq(n, d=dt)
    domega = (2*np.pi)/(dt*(n + 1))
    omegas = 2 * np.pi * freq
    omegas = omegas[1:]

    # Fast Fourier Transform of Acceleration
    accfft = np.array(np.fft.fft(acc, axis=0))

    # High-pass filter
    n_low = int(np.floor((2*np.pi*f_low)/domega))
    accfft[:n_low - 1] = 0
    accfft[n - n_low:] = 0

    # Integrating the Fast Fourier Transform
    velfft = -accfft[1:] / (omegas * 1j)
    posfft = accfft[1:] / ((omegas * 1j) ** 2)

    velfft = np.array([0] + list(velfft))
    posfft = np.array([0] + list(posfft))

    # Inverse Fast Fourier Transform back to time domain
    vel = np.real(np.fft.ifft(velfft, axis=0))
    pos = np.real(np.fft.ifft(posfft, axis=0))

    return vel, pos

This answer was posted as an edit to the question Different results using FFT in Matlab compared to Python by the OP Kakemonster under CC BY-SA 3.0.

Improve this answer

1 Comment

Consider making this a community wiki.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.