-1

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')
4

2 回答 2

1

I am adding another answer in reply to your followup questions. A rectangular window (or boxcar window) of width T, when used as moving average filter, is a lowpass filter with transfer function magnitude

|H(f)|=(sin(pifT))/(pifT) (1)

as you can verify by doing the Fourier transform. Therefore H(f)=0.707 when fco=0.44/T, and H(f)=0 when f=1/T. So a lowpass filter with cutoff fco can be implemented with a rectangular window with width on the order of 1/fco (specifically .44/fco). When fco=0.03 Hz, width is on the order of 30 s (specifically 15 s). See plot of transfer function of a 30 s-wide window. Transfer function of a rectangular window of width 30 s.

How to find the necessary order:

>> Fs=200; Fc=.03;
>> Hd003=designfilt('lowpassfir','CutoffFrequency',Fc,'SampleRate',Fs);
Hd003 = designfilt('lowpassfir', 'PassbandFrequency', .03, 'StopbandFrequency', .1, 'PassbandRipple', 3, 'StopbandAttenuation', 20, 'SampleRate', 200, 'DesignMethod', 'kaiserwin');
>> disp(length(Hd003.Coefficients));
        2400

The command Hd003=designfilt(...) above causes the Filter Design Assistant to launch, because designfilt() does not have enough info. In the FDA, specify Order mode = Minimum, Passband freq=0.03, Stopband freq=0.1, passband ripple=3 dB, stopband atten.=20 db, design method=Kaiser window. I chose those frequencies and dB values because a 2nd order Butterworth does equally well. (See figure.) Then click "OK", and control returns to the Matlab command window, which runs 'designfilt' with the specified parameters. Then check the length of the coefficient vector: it is 2400! I.e. the FIR filter order=2399, and the filter is 12 seconds long. This is not practical, and it is totally un-usable on signals as short as your 1.6 second long example. Filter Design Assistant screenshot Freq. response of a 2nd order Butterworth lowpass with Fc=0.03 Hz. From Matlab Filter Designer. Freq. response of FIR filter meeting design criteria similar to a 2nd order Butterworth.  The order required was 2399.  From Matlab Filter Designer.

  1. This equation for |H(f)| is for continuous time. It is accurate for discrete time also, as long as the sampling rate Fs>=10/T.
于 2021-01-18T03:31:54.197 回答
0

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.

于 2021-01-17T09:03:10.287 回答