0

我对加速度数据使用了巴特沃斯带通滤波器。现在我想对齐/组合原始信号和平滑信号以进行观察/分析。问题是我正在尝试对齐，但它不起作用。可能我做错了。我需要指导，如何调整它们。下面我发布了我试图用来对齐两个信号的输出图像和代码。我放在这里的输出图像由原始信号和使用巴特沃斯带通滤波器完成的滤波信号组成。谢谢

delay = mean(grpdelay(B,A))
tt = Time(1:end-delay);
sn = norm(1:end-delay);

sf = Data_filtered;
sf(1:delay) = [];
plot(tt,sn)
hold on, plot(tt,sf,'-r','linewidth',1.5), hold off
title 'Gyro Filtration'
xlabel('Time (s)'), legend('Original Signal','Filtered Shifted Signal')

输出原始信号和滤波信号

4

1 回答 1

0

我建议您使用 filtfilt() 进行过滤，因为 filfilt() 在所有频率上都具有完全零延迟。它实现了这一点，因为它向前过滤信号，然后向后过滤。向前通过时的时间延迟被向后通过时相等但相反的延迟抵消。请注意，由于您进行了两次过滤，因此每个频率的最终衰减（如果以 dB 为单位）是您从单次通过中获得的衰减量的两倍。因此，您会在原始截止频率处获得 6 dB 的衰减，而不是通常的 3 dB。对于大多数应用程序，包括您的应用程序，我相信这不是问题。如果这是一个问题，那么您可以调整初始滤波器定义的截止值，以便在两次通过后在所需截止值处获得 3 dB 衰减。

有关延迟与 filt() 和 filtfilt() 的比较，请参见随附的代码和图。该图显示 raw 和 filtfilt() 信号都在完全相同的时间通过 y=0 转换，这反映了 filtfilt() 没有任何延迟。该图还显示，原始信号中的高频噪声已被滤波器去除。

比较原始信号和 filter() 和 filtfilt() 的输出，对于 Fhpc=0.1 Hz、Flpc=5 Hz 的 4 阶带通 Butterworth。

如果由于某种原因您不想使用 filtfilt()，那么通带中 n 阶低通巴特沃斯滤波器的延迟的一个非常好的近似值是

延迟 = Kn/Fc（延迟单位为 s，Fc=截止频率单位为 Hz）

其中 Kn 是取决于滤波器阶数 n 的常数：

n=2,3,4,6,8

Kn=0.225,0.318,0.416,0.615,0.816

有关更多命令和推导，请参阅我的论文“A general solution for the time delay...”，J Biomechanics (2007)，https://pubmed.ncbi.nlm.nih.gov/16545388/。

%filterCompare.m  WCR 20210118
% Question on stack exchange
clear;

%Make an illustrative input signal
Fs=200;                     %sampling frequency (Hz)
N=800;                      %number of points
Tdur=N/Fs;                  %duration (s)
Time=(0:(N-1))/Fs;          %time vector
x=square(Time*2*pi/Tdur)...     %square wave plus high frequency ripple
    +0.2*sin(2*pi*.25*Fs*Time);

%Compute the filter coefficients
order=4;
Fhpc=0.1;                   %high pass cutoff frequency (Hz)
Flpc=5;                     %low pass cutoff frequency (Hz)
Wn1=[Fhpc,Flpc]*2/Fs;       %normalized cutoff frequencies
[b,a]=butter(order,Wn1,'bandpass');  %compute filter coefficients

%Filter the data
y1=filter(b,a,x);         %standard (forward-only) filter
y2=filtfilt(b,a,x);       %forward and backward filter

%Plot the data
subplot(2,1,1);
plot(Time,x,'r.-',Time,y1,'g.-',Time,y2,'b.-');
ylabel('Amplitude');
legend('Raw Data','filter()','filtfilt()');
grid on;
subplot(2,1,2);
plot(Time,x,'r.-',Time,y1,'g.-',Time,y2,'b.-');
xlabel('Time (s)'); ylabel('Amplitude');
legend('Raw Data','filter()','filtfilt()');
axis([.45*Tdur .55*Tdur -inf inf]);
grid on;
于 2021-01-18T20:01:33.043 回答