Main Content

filtfilt

零相位数字滤波

说明

y = filtfilt(b,a,x) 通过正向和反向处理输入数据 x 来执行零相位数字滤波。在对数据进行正向滤波后,该函数会匹配初始条件以最大限度地减少启动和结束时的瞬变,反转滤波后的序列,并让反转后的序列再次通过滤波器。结果具有以下特征:

  • 零相位失真

  • 滤波器传递函数等于原始滤波器传递函数的平方幅值

  • 滤波器阶数是由 ba 指定的滤波器阶数的两倍

filtfilt 实现了 Gustafsson 提出的算法 [1]

不要将 filtfilt 与微分器和希尔伯特 FIR 滤波器结合使用,因为这些滤波器的运算在很大程度上取决于它们的相位响应。

y = filtfilt(sos,g,x) 使用由矩阵 sos 和定标值 g 表示的二阶节 (biquad) 滤波器对输入数据 x 进行零相位滤波。

示例

y = filtfilt(d,x) 使用数字滤波器 d 对输入数据 x 进行零相位滤波。使用 designfilt 根据频率响应设定生成 d

示例

全部折叠

零相位滤波有助于将滤波后的时间波形中的特征准确地保留在它们在未经滤波的信号中出现的位置。

对合成心电图 (ECG) 波形进行零相位滤波。生成波形的函数在示例的末尾。QRS 复波是 ECG 的一个重要特征。如下图所示,它在时间点 160 左右开始。

wform = ecg(500);

plot(wform)
axis([0 500 -1.25 1.25])
text(155,-0.4,"Q")
text(180,1.1,"R")
text(205,-1,"S")

用加性噪声损坏 ECG。重置随机数生成器以获得可重现的结果。构造一个低通 FIR 等波纹滤波器,并使用零相位滤波和传统滤波对含噪波形进行滤波。

rng default

x = wform' + 0.25*randn(500,1);
d = designfilt("lowpassfir", ...
    PassbandFrequency=0.15,StopbandFrequency=0.2, ...
    PassbandRipple=1,StopbandAttenuation=60, ...
    DesignMethod="equiripple");
y = filtfilt(d,x);
y1 = filter(d,x);

subplot(2,1,1)
plot([y y1])
title("Filtered Waveforms")
legend("Zero-phase Filtering","Conventional Filtering")

subplot(2,1,2)
plot(wform)
title("Original Waveform")

零相位滤波可减少信号中的噪声,并在原始信号中出现 QRS 复波的同时保留该复波。传统滤波可减少信号中的噪声,但会使 QRS 复波延迟。

使用巴特沃斯二阶节滤波器重复进行滤波。

d1 = designfilt("lowpassiir",FilterOrder=12, ...
    HalfPowerFrequency=0.15,DesignMethod="butter");
y = filtfilt(d1,x);

subplot(1,1,1)
plot(x)
hold on
plot(y,LineWidth=3)
legend("Noisy ECG","Zero-Phase Filtering")

此函数生成 ECG 波形。

function x = ecg(L)
%ECG Electrocardiogram (ECG) signal generator.
%   ECG(L) generates a piecewise linear ECG signal of length L.
%
%   EXAMPLE:
%   x = ecg(500).';
%   y = sgolayfilt(x,0,3); % Typical values are: d=0 and F=3,5,9, etc. 
%   y5 = sgolayfilt(x,0,5); 
%   y15 = sgolayfilt(x,0,15); 
%   plot(1:length(x),[x y y5 y15]);

%   Copyright 1988-2002 The MathWorks, Inc.

a0 = [0,1,40,1,0,-34,118,-99,0,2,21,2,0,0,0]; % Template
d0 = [0,27,59,91,131,141,163,185,195,275,307,339,357,390,440];
a = a0/max(a0);
d = round(d0*L/d0(15)); % Scale them to fit in length L
d(15)=L;

for i=1:14
       m = d(i):d(i+1)-1;
       slope = (a(i+1)-a(i))/(d(i+1)-d(i));
       x(m+1) = a(i)+slope*(m-d(i));
end

end

输入参数

全部折叠

传递函数系数,指定为向量。如果使用全极点滤波器,请为 b 输入 1。如果使用全零 (FIR) 滤波器,请为 a 输入 1

示例: b = [1 3 3 1]/6a = [3 0 1 0]/3 用于指定一个归一化 3 dB 频率为 0.5π 弧度/采样点的三阶巴特沃斯滤波器。

数据类型: single | double

输入信号,指定为实数值或复数值向量、矩阵或 N 维数组。x 包含的值必须为有限值。x 的长度必须大于滤波器阶数的三倍,定义为 max(length(B)-1, length(A)-1)。除非 x 是行向量,否则该函数沿 x 的第一个数组维度进行运算。如果 x 是行向量,则该函数沿第二个维度进行运算。

示例: cos(pi/4*(0:159))+randn(1,160) 是单通道行向量信号。

示例: cos(pi./[4;2]*(0:159))'+randn(160,2) 是双通道信号。

数据类型: single | double
复数支持:

二阶节系数,指定为矩阵。sos 是一个 K×6 矩阵,其中节数 K 必须大于或等于 2。如果节数小于 2,则该函数将输入视为分子向量。sos 的每行对应于二阶 (biquad) 滤波器的系数。sos 的第 i 行对应于 [bi(1) bi(2) bi(3) ai(1) ai(2) ai(3)]

示例: s = [2 4 2 6 0 2;3 3 0 6 0 0] 用于指定一个归一化 3 dB 频率为 0.5π 弧度/采样点的三阶巴特沃斯滤波器。

数据类型: single | double

缩放因子,指定为向量。

数据类型: single | double

数字滤波器,指定为 digitalFilter 对象。使用 designfilt 根据频率响应设定生成数字滤波器。

示例: d = designfilt("lowpassiir",FilterOrder=3,HalfPowerFrequency=0.5) 用于指定一个归一化 3 dB 频率为 0.5π 弧度/采样点的三阶巴特沃斯滤波器。

输出参量

全部折叠

滤波后的信号,以向量、矩阵或 N 维数组形式返回。

如果 filtfilt 的输入是单精度值,则该函数返回单精度的输出 y

参考

[1] Gustafsson, F. “Determining the initial states in forward-backward filtering.” IEEE® Transactions on Signal Processing. Vol. 44, April 1996, pp. 988–992. https://doi.org/10.1109/78.492552.

[2] Mitra, Sanjit K. Digital Signal Processing. 2nd Ed. New York: McGraw-Hill, 2001.

[3] Oppenheim, Alan V., and Ronald W. Schafer, with John R. Buck. Discrete-Time Signal Processing. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.

扩展功能

版本历史记录

在 R2006a 之前推出