Designing Filter for Accelerometer

Hello everyone,

I want to design filter for my accelerometer to get position. The objective is to get the trace in from the X and Y axes. I have tried to code some in matlab but i still cant get it. I have attached the excel data and matlab code that I have prepared. The excel data is a diagonal of approximately 66cm with 50cm and 44cm for X and Y respectively. I really need help for this because this is my first time designing filter for very noisy sensor. Any filter will do.

Thank you in advance!

Answer 1

[Star Strider]

The data are not regularly-sampled, so in order to do any meaningful signal processing on them, they need to be resampled to a consistent sampling interval. I did this with the resample function, and since the sampling frequency of the original data was about 73 Hz, I chose a sampling frequency of 70 Hz. Choose whatever works best, however it would be better to go with a lower than a higher resampling frequency, since this avoids creating data where none previously existed.

The noise ais actually broadband, and the best way to deal with broadband noise is to use the Savitzky-Golay filter (sgolayfilt function) or wavelet denoising. I will help you with sgolayfilt if that is the direction you want to go. I use a frequency selective filter here to start with.

% opts = weboptions('ContentType','text');
% W = webread('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1038845/Tracking.m', opts)

T1 = readtable('https://www.mathworks.com/matlabcentral/answers/uploaded_files/1038850/Diagonal66cm.xlsx', 'VariableNamingRule','preserve')

T1 = 1696×4 table
     AcX(G)      AcY(G)     AcZ(G)    Timestamps
    ________    ________    ______    __________

     0.05481    0.019714    1.0432          0   
    0.060181    0.016541    1.0424      0.013   
    0.056519    0.014343    1.0468      0.024   
    0.058716    0.021179     1.041      0.039   
    0.053833    0.013855    1.0344       0.05   
    0.055298    0.015076    1.0439      0.061   
    0.060181    0.015076    1.0495      0.072   
    0.060669    0.025818    1.0505      0.083   
    0.053589    0.021667    1.0539      0.094   
    0.059448    0.016785     1.048      0.105   
    0.056763    0.020447    1.0424      0.116   
    0.061157    0.012878    1.0434      0.131   
    0.049927    0.017029    1.0475      0.142   
    0.049194    0.013855    1.0446      0.153   
    0.052368    0.021179    1.0485      0.164   
    0.059692    0.018005    1.0458      0.175   

tstr = {'X','Y','Z'};

figure
for k = 1:3
    subplot(3,1,k)
    plot(T1.Timestamps, T1{:,k})
    grid
    title(tstr{k})
    xlim([min(T1.Timestamps) max(T1.Timestamps)])
end
xlabel('Time')
sgtitle('Original Accelerometer Data')

fse = 1/mean(diff(T1.Timestamps))

fse = 73.5869

Fs = 70;
[Acr,tr] = resample(T1{:,1:3}, T1.Timestamps, Fs)

Acr = 1613×3

0.0538 0.0193 1.0185 0.0603 0.0164 1.0520 0.0568 0.0162 1.0408 0.0578 0.0199 1.0418 0.0541 0.0130 1.0382 0.0604 0.0179 1.0500 0.0583 0.0246 1.0510 0.0550 0.0195 1.0527 0.0594 0.0177 1.0434 0.0578 0.0165 1.0450

tr = 1613×1

0 0.0143 0.0286 0.0429 0.0571 0.0714 0.0857 0.1000 0.1143 0.1286

figure
for k = 1:3
    subplot(3,1,k)
    plot(tr, Acr(:,k))
    grid
    title(tstr{k})
    xlim([min(tr) max(tr)])
end
xlabel('Time')
sgtitle('Resampled Accelerometer Data')

Fn = Fs/2;
L = numel(tr);

NFFT = 2^nextpow2(L);                                               % For Efficiency
FT_Acr = fft(Acr-mean(Acr),NFFT)/L;                                 % Subtract 'mean' TO See Other Peaks
Fv = linspace(0, 1, NFFT/2+1)*Fn;                                   % Frequency Vector
Iv = 1:numel(Fv);                                                   % Index Vector

figure
for k = 1:3
    subplot(3,1,k)
    plot(Fv, abs(FT_Acr(Iv,k))*2)
    grid
    title(tstr{k})
    xlim([min(Fv) max(Fv)])
end
xlabel('Frequency')
ylabel('Amplitude')
sgtitle('Accelerometer: Fourier Transform')
xlim([min(Fv) max(Fv)])

% return

Wp = [0.00005 12.5]/Fn;                                             % Passband Frequency (Normalised)
Ws = [0.9 1.10].*Wp;                                                % Stopband Frequency (Normalised)
Rp =  1;                                                            % Passband Ripple
Rs = 60;                                                            % Passband Ripple (Attenuation)
[n,Wp] = ellipord(Wp,Ws,Rp,Rs);                                     % Elliptic Order Calculation
[z,p,k] = ellip(n,Rp,Rs,Wp);                                        % Elliptic Filter Design: Zero-Pole-Gain 
[sos,g] = zp2sos(z,p,k);                                            % Second-Order Section For Stability

figure
freqz(sos, 2^16, Fs)                                                % Filter Bode Plot

Acr_filt = filtfilt(sos, g, Acr);                                   % Filter Signal
tstr = {'X','Y','Z'};

figure
for k = 1:3
    subplot(3,1,k)
    plot(tr, Acr_filt(:,1))
    grid
    title(tstr{k})
    xlim([min(tr) max(tr)])
end
xlabel('Time')
sgtitle('Filtered Accelerometer Data')

.