为什么即使在更改加速时间序列数据的截止频率后,fir 滤波器也能得到相同的输出?

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

我对过滤有疑问。

目前我正在使用惯性测量单元(accelerometer/acceleration数据)进行行人步数检测。我需要在预处理级别过滤我的信号。谁能建议哪一个是获得良好准确性的最佳过滤算法? 现在我使用 digital lowpass fir filter with order=2, window size=1.2sec,sampling frequecy=200hz 但似乎没有在职的。我想使用较低的截止频率。我使用了 0.03hz 和 3hz 的截止频率 ,但我得到了 两个截止频率相同的结果 。我需要您的指导或帮助,我该如何进行。下面,我分别在 3hz 和 0.03hz 附上过滤结果的图像 link 以及我尝试过的代码。 有人可以建议我或在 matlab 中提供任何好的过滤器 and/or 我该如何处理?

fir result at 3hz cutoff fir result at 0.03hz cutoff

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

您试图制作二阶 FIR 滤波器。对于您的采样率 (200 Hz) 和所需的截止频率(3 或 0.03 Hz),该阶数太低了。 FIR 滤波器所需的阶数与 IIR 滤波器的阶数有很大不同。 N 阶 FIR 滤波器对 N+1 个数据点进行移动平均。 Hd=designfilt(...) 计算移动平均值的系数或权重。我使用您的代码片段制作了 3 Hz 和 0.03 Hz 截止滤波器,并将它们命名为 Hd300 和 Hd003。然后我查看了两个滤波器的系数:

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

如您所见,它们实际上是相同的——这就是输出滤波信号看起来相同的原因。系数非常相似,因为当采样率为 200 Hz 时,order=2 太低而无法制作 3 或 0.03 Hz 截止频率的有效 FIR 滤波器。 您尝试的 0.03 截止频率对应于大约 30 (=1/.03) 秒的时间常数。对只有 1.6 秒长的数据使用如此低的截止值是没有意义的。即使你有数百秒的数据,它也没有意义,因为你会用大约 30 秒的宽度 window 来平滑数据,这将很难检测到每一步。 更好的方法:简单的二阶巴特沃斯低通,截止频率 = 1 到 10 Hz。请参见下面的代码并参见图。我在我的代码中使用了 filtfilt() 而不是 filter(),因为 filtfilt() 处理启动瞬态更好。将 filtfilt() 替换为 filter() 你就会明白我的意思了。

%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.

我正在添加另一个答案来回复您的后续问题。宽度为 T 的矩形 window(或 boxcar window)用作移动平均滤波器时,是具有传递函数幅度

的低通滤波器

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

你可以通过傅里叶变换来验证。因此当fco=0.44/T时H(f)=0.707,当f=1/T时H(f)=0。因此,可以使用宽度约为 1/fco(具体为 .44/fco)的矩形 window 来实现具有截止 fco 的低通滤波器。当 fco=0.03 Hz 时,宽度约为 30 s(具体为 15 s)。请参见 30 秒宽的传递函数图 window。

如何找到需要的订单:

>> 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

上面的命令 Hd003=designfilt(...) 会启动滤波器设计助手,因为 designfilt() 没有足够的信息。在 FDA 中,指定 Order mode = Minimum,Passband freq=0.03,Stopband freq=0.1,passband ripple=3 dB,stopband atten.=20 db,design method=Kaiser window。我之所以选择这些频率和 dB 值,是因为二阶 Butterworth 的表现同样出色。 (见图)然后点击“确定”,控制returns到Matlab命令window,按指定参数运行'designfilt'。然后查看系数向量的长度:是2400! IE。 FIR 滤波器阶数=2399,滤波器长度为 12 秒。这是不切实际的,而且它在短至 1.6 秒长示例的信号上完全不可用。

  1. 此等式为 |H(f)|是连续时间。对于离散时间也是准确的,只要采样率 Fs>=10/T.