Search code examples
matlabmachine-learningfilteringaccelerometerbutterworth

Why getting same output from fir filter even after changing the cutoff frequency for acceleration time series data?


I have question about filtering.

Currently I am working on pedestrian step detection using an inertial measuring unit (accelerometer/acceleration data). I need to filter my signal at preprocessing level. Could anyone suggest which one will be the best filtration algorithm to get good accuracy? For now I used digital lowpass fir filter with order=2, window size=1.2sec,sampling frequecy=200hz but seems not working. I want to use a lower cutoff frequency. I used 0.03hz and 3hz of cutoff frequecny but I am getting the same result for both cutoff frequencies. I need your guidance or help that how can I proceed. Below, I am attaching the images of filtration result as link at 3hz and 0.03hz respectively and also code that I tried. Could some one suggest me or provide any good filter in matlab and/or how can I work on that?

fir result at 3hz cutoff fir result at 0.03hz cutoffFull signal

Fs = 200;
Hd = designfilt('lowpassfir','FilterOrder',2,'CutoffFrequency',0.03, ...
   'DesignMethod','window','Window',{@kaiser,1.2},'SampleRate',Fs);
y1 = filter(Hd,norm);

plot(Time,norm,'b', Time, y1,'red')

xlabel('Time (s)')
ylabel('Amplitude')
legend('norm','Filtered Data')

Solution

  • You tried to make a 2nd order FIR filter. That order is far to low for your sampling rate (200 Hz) and desired cutoff frequency (3 or 0.03 Hz). The required order in a FIR filter is very different from the order of an IIR filter. An Nth order FIR filter does a moving average of N+1 data points. Hd=designfilt(...) computes the coefficients, or weights, for the moving average. I made the 3 Hz and 0.03 Hz cutoff filters, using your code fragment, and called them Hd300 and Hd003. Then I viewed the coefficients for the two filters:

    >> disp(Hd003.Coefficients);
        0.2947    0.4107    0.2947
    >> disp(Hd300.Coefficients);
        0.2945    0.4110    0.2945
    

    As you can see, they are practically identical - which is why the output filtered signals look the same. The coefficients are very similar because order=2 is far, far too low to make an effective FIR filter of 3 or 0.03 Hz cutoff, when the sampling rate is 200 Hz. The 0.03 cutoff frequency which you tried corresponds to a time constant of about 30 (=1/.03) seconds. It does not make sense to use such a low cutoff on data which is only 1.6 seconds long. Even if you had hundreds of seconds of data, it would not make sense, because you would be smoothing the data with about a 30 second wide window, which would make it very hard to detect each step. A better approach: a simple 2nd order Butterworth lowpass, cutoff frequency=1 to 10 Hz. See code below and see figure. I used filtfilt() in my code rather than filter(), because filtfilt() handles the startup transient better. Replace filtfilt() with filter() and you will see what I mean.

    %pedometerAccelFilter.m  WCR 20210117
    % Question on stack exchange
    % The code provided on stack exchange assumes vector "norm"
    % contains 1.6 seconds of acceleration data sampled at 200 Hz.
    % I digitized the data from the figure on stack exchange and
    % saved it in file pedometerInterpData.txt (2 columns, 329 rows).
    
    %Load the data
    data=load('pedometerInterpData.txt');
    Time=data(:,1); norm=data(:,2);
    
    %Compute the filter coefficients
    Fs=200;                     %sampling frequency (Hz)
    order=2;
    Fc1=5;                      %filter 1 cutoff frequency (Hz)
    Wn1=Fc1*2/Fs;               %filter 1 normalized cutoff frequency
    [b1,a1]=butter(order,Wn1);  %filter 1 coefficients
    Fc2=2;                      %filter 2 cutoff frequency (Hz)
    Wn2=Fc2*2/Fs;               %filter 2 normalized cutoff frequency
    [b2,a2]=butter(order,Wn2);  %filter 2 coefficients
    
    %Filter the data
    y1=filtfilt(b1,a1,norm);    %filtfilt() could be replaced with filter()
    y2=filtfilt(b2,a2,norm);
    
    %Plot the data
    plot(Time,norm,'r.-',Time,y1,'gx-',Time,y2,'bo-');
    xlabel('Time (s)'); ylabel('Amplitude');
    legend('Raw Data','Filter 1','Filter 2');
    

    Figure created by code above.