Signal Processing For Integration

Hello, i'm doing an experimental Program. I read some data from an IMU sensor, precisely i read accelerometer data along 3 axis (x, y, z), every 0,04s.

Now i want to pass from vertical acceleration data (in my case the vertical acceleration is on the y-axis), to vertical displacement data. I know the noisy of signal, i read lot about different double integration technique. I tried Euler Method, but the drift is very high (and I know that this is a consequence of using this method).

I read instead of the Kalman filter, or by using Fourier offer approximations more real.

I want to focus particularly on Fourier, are not a numerical analysis expert, i wanted to ask if someone can give me some advice on how to set it up, and then switch between the time domain to a frequency domain, and then do the double integration to obtain displacement and return in the time domain.

Can someone help me?

Answer 1

[Star Strider]

The documentation for the fast Fourier transform is here: fft (link). Particularly note the code between the first (top) two plot figures. Use the Fourier transform to determine what part of the spectrum are your signals, and what part are noise.

First, I would use a bandpass filter to eliminate the constant (d-c) offset at the low end and the noise at the high end. In my opinion, FIR filters are best for this. The easiest way to design your filter is to use the designfilt function.

Always check the performance of your filter with the freqz function to be certain it is doing what you want it to do.

Use the filtfilt function to do the actual filtering. It has a maximally-flat phase characteristic, so no phase distortion will appear in your filter regardless of the filter design you use. (You do not have to use a Bessel filter for this.)

Second, after you have filtered your signals, integrate them with the cumtrapz (or trapz) function. The reason to do the filtering first is to remove the constant offset of low-frequency baseline variation, because if you do not, these will integrate with your signal, producing meaningless results for your displacement data.

[sbareben]

Thanks for the answer @Star Strider, i tried to use fft and pay particullary attention to code between the first two plo figures as you suggest.

I change the time period of the register to 10ms.

I post a sample of my original registered data in the fig (registeredAcc.fig).

After i read the original data, i've apply a code to remove gravity.

function[AccelerometerY_WithouthGravity] =removeGravity(RawAccelerometerY, Time)
%https://developer.android.com/reference/android/hardware/SensorEvent.html#values in merita
num_elements = numel(Time);
alpha = zeros(num_elements,1);
for index = 1:numel(Time)-1
  alpha(index) = Time(index)/(Time(index) + (Time(index+1)-Time(index)));
end
alpha(numel(alpha), 1 ) = Time(numel(alpha))/(Time(numel(alpha)) + (Time(numel(alpha))-Time(numel(alpha)-1)));
gravity = zeros(num_elements,1);
gravity(1,1) = alpha(1)*0 + (1-alpha(1))*RawAccelerometerY(1);
for index = 2:numel(RawAccelerometerY)
   gravity(index, 1) = alpha(index)*gravity(index-1) + (1-alpha(index))*RawAccelerometerY(index);
end
AccelerometerY_WithouthGravity = zeros(num_elements,1);
for index = 1:numel(gravity)
  AccelerometerY_WithouthGravity(index,1) =  RawAccelerometerY(index) - gravity(index);
end

i've used this beacuse i read from a smartphone and i've watch from android developer how to remove it. (i hope this is correct)

And the result i obtain i show in the fig (NoGravity.fig)

Than i tried to use fft like this:

Fs =0.1;            % Sampling frequency KHZ                   
T = 1/Fs;           % Sampling period ms        
L = 2000;             % Length of signal 
t = (0:L-1)*T;        % Time vector 
Y = fft(AccelerometroY);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
plot(f,P1)

The result i obtain is in the figure (fft.fig).

Now if this is correct (but i'm not really sure), to find amplitude (precisly i'm interested to high peak of signal) i need to filter my data, to do it, i need to use the filter you suggest me?

[Star Strider]

My pleasure.

Your Fourier analysis code appears correct to me. Instead of plotting only ‘P1(2:end-1)’ (this shifts the signal with respect to the frequency axis), subtract the mean of your signal before you take the fft, then plot all of it.

For example:

Y = fft(AccelerometroY - mean(AccelerometroY));
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
f = Fs*(0:(L/2))/L;
plot(f,P1)

I did not test this. It should work.

This filter is stable, and I believe does what you want. Change only the higher passband and stopband frequencies if you want to change its properties. (The freqz call shows the filter properties. It is not necessary for the rest of the code.)

The Code —

Fs = 0.1;
Fn = Fs/2;                                              % Nyquist Frequency
Wp = [2E-4  3E-3]/Fn;                                   % Passband Frequencies (Normalized)
Ws = [1E-4  4E-3]/Fn;                                   % Stopband Frequencies (Normalized)
Rp = 10;                                                % Passband Ripple (dB)
Rs = 50;                                                % Stopband Ripple (dB)
[n,Ws] = cheb2ord(Wp,Ws,Rp,Rs);                         % Filter Order
[c,b,a] = cheby2(n,Rs,Ws);                              % Filter Design
[sosbp,gbp] = zp2sos(c,b,a);                            % Convert To Second-Order-Section For Stability

figure(1)
freqz(sosbp, 2^16, Fs)                                  % Bode Plot Of Filter
set(subplot(2,1,1), 'XLim',[0  5E-3])                   % ‘Zoom’ X-Axis To See Passband
set(subplot(2,1,2), 'XLim',[0  5E-3])                   % ‘Zoom’ X-Axis To See Passband

Then filter your signal with the filtfilt function:

AccelerometroY_filtered = filtfilt(sosbp, gbp, AccelerometroY);

[sbareben]

Thanks again, i test the filter and if i changing the passband and stopband frequencies like i need he should work well, but i would ask you just another think, than i filter the signal i integrate twice:

vel = cumtrapz(AccelerationY_Filtered);
disp = cumtrapz(vel);

The result of integration is plotted on the fig.

I do not understand why the displacement to be increasing, and doesn't assume a pattern as uniform as possible on the horizontal axis, I am interested only to high peaks of acceleration and I expect to see the same thing at vertical displacement, and after integration seems to me that the less coincide.

[Star Strider]

My pleasure.

The displacement is increasing (rather than oscillating) probably because you are integrating a constant or low-frequency baseline variation. You may have to increase the lower passband and stopband frequencies in your filter to elinimiate a low-frequency trend. (I did linear regression detrending on the filtered acceleration signal and got an improved result.)

You may have to experiment with the passband and stopband frequencies to get the result you want.

[sbareben]

ok, i try to experiment better, to integrate the signal the function cumtrapz is the best i can do? i read about omega algo, do you suggest me anything else?

[Star Strider]

The cumtrapz function is likely your only option. I do not recognize ‘omega algo’, so I cannot comment on it.

You appear to be doing everything correctly.

[sbareben]

Thank you again, if I wanted to get exclusively the high peaks of my signal, and eliminate all that lies below a frequency, so let's get a signal with low peaks tending say to 0, and the high unaffected peaks and then only rounded, should I change the filter?

[Star Strider]

My pleasure.

I am not certain that I understand what you want to do.

Your best option may be the findpeaks function if you want to choose specific peaks by height or other criteria.

[sbareben]

Sorry again, but in previous sample code of Cheby filter that you post, how did you find or set the stopband and passband ripple? Are calculated in function of Passband and stopband frequencies? And how do you determine the Passband and Stopband frequencies?

[Star Strider]

The passband and stopband ripple amplitudes are essentially just ‘what works’. In a Chebyshev Type II filter, the passband ripple is essentially irrelevant, because it has a flat passband. A large value usually results in a stable, short filter. The stopband ripple is the attenuation in the stopband, so a good choice is 50 dB, with a stopband amplitude of 1E-5.

Filter design is a compromise between the performance the filter designer wants, and the practical considerations in the filter implementation. This applies to both continuous and discrete filters. In continuous filters, the problem is hardware complexity, in discrete filters, the problem is filter length.