如何对齐原始信号和过滤后的信号进行分析?

How to align original and filtered signal for analysis?

我使用 巴特沃斯带通滤波器 作为我的加速度数据。现在我想 align/combine 原始信号和平滑信号 用于 obervation/analysis . 问题是我试图对齐但它不起作用。可能我做错了。我需要指导,如何才能对齐它们。 下面我将发布我的输出图像和代码,我试图使用它来 align 这两个信号。我放在这里的输出图像由 原始信号和使用 巴特沃斯带通滤波器 完成的滤波信号 组成。谢谢

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

output original and filtered signals

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

filt() 和 filtfilt() 的延迟比较见附件代码和图。该图显示 raw 和 filtfilt() 信号同时通过 y=0 转换,这反映了 filtfilt() 没有任何延迟。该图还表明,原始信号中的高频噪声已被滤波器去除。

如果出于某种原因您不想使用 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

有关更多命令和推导,请参阅我的论文“时间延迟的一般解决方案...”,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;