filtfilt

Filtrage numérique à phase zéro

Syntaxe

y = filtfilt(b,a,x)

y = filtfilt(sos,g,x)

y = filtfilt(d,x)

y = filtfilt(B,A,x,"ctf")

y = filtfilt({B,A,g},x)

Description

y = filtfilt(b,a,x) applique un filtrage numérique à phase zéro en traitant les données d’entrée x à la fois en sens avant et arrière. Après avoir filtré les données en sens avant, la fonction reproduit les conditions initiales pour minimiser les transitoires de démarrage et de fin, inverse la séquence filtrée et refait passer la séquence inversée par le filtre. Le résultat présente les caractéristiques suivantes :

• Distorsion de phase nulle

• Fonction de transfert de filtre égale à l’amplitude au carré de la fonction de transfert de filtre d’origine

• Ordre de filtre égal au double de l’ordre du filtre défini par b et a

filtfilt implémente l’algorithme proposé par Gustafsson [1].

N’utilisez pas filtfilt avec des filtres FIR différenciateurs et de Hilbert, car le fonctionnement de ces filtres dépend fortement de leur réponse en phase.

y = filtfilt(sos,g,x) applique un filtrage à phase zéro aux données d’entrée x, en utilisant le filtre à sections de second ordre (biquadratique) représenté par la matrice sos et les valeurs d’échelle g.

y = filtfilt(d,x) applique un filtrage à phase zéro aux données d’entrée x, en utilisant un filtre numérique d. Utilisez designfilt pour générer d selon des spécifications de réponse en fréquence.

y = filtfilt(B,A,x,"ctf") applique un filtrage à phase zéro aux données d’entrée x, en utilisant des Cascaded Transfer Functions (CTF) définies par les coefficients de numérateur B et de dénominateur A. (depuis R2024b)

Remarque

Spécifiez l’option "ctf" pour distinguer les matrices de numérateurs CTF B à six colonnes des matrices de sections de second ordre en entrée sos lorsque vous définissez A sous forme de scalaire ou de vecteur.

y = filtfilt({B,A,g},x) inclut les valeurs scalaires g du filtre dans le filtrage des données d’entrée x. (depuis R2024b)

Exemples

Filtrage à phase zéro d’une forme d’onde d’électrocardiogramme

Le filtrage à phase zéro permet de conserver les caractéristiques d’une forme d’onde temporelle filtrée en les faisant apparaître exactement au même endroit que dans le signal non filtré.

Appliquez un filtre à phase zéro à une forme d’onde d’électrocardiogramme (ECG) synthétique. La fonction qui génère la forme d’onde se trouve à la fin de l’exemple. Le complexe QRS est une caractéristique importante de l’ECG. Ici, il commence vers le point temporel 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")

Altérez l’ECG avec du bruit additif. Réinitialisez le générateur de nombres aléatoires pour obtenir des résultats reproductibles. Créez un filtre FIR passe-bas à ondulations constantes et filtrez la forme d’onde bruitée en utilisant à la fois un filtrage à phase zéro et conventionnel.

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

tiledlayout("flow")
nexttile
plot([y y1])
title("Filtered Waveforms")
legend(["Zero-phase Filtering" "Conventional Filtering"])

nexttile
plot(wform)
title("Original Waveform")

Le filtrage à phase zéro réduit le bruit dans le signal et fait apparaître le complexe QRS au même moment que dans le signal d’origine. Le filtrage conventionnel réduit le bruit dans le signal mais retarde le complexe QRS.

Répétez le filtrage, cette fois avec un filtre de Butterworth à sections de second ordre.

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

figure
plot(x)
hold on
plot(y,LineWidth=3)
hold off
legend(["Noisy ECG" "Zero-Phase Filtering"])

La fonction suivante génère la forme d’onde de l’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

Filtrage en cascade à phase zéro d’un chirp linéaire

Depuis R2024b

Utilisez des fonctions de transfert en cascade pour appliquer un filtrage à phase zéro.

Créez un filtre elliptique avec une ondulation dans la bande passante de 0,1 dB et une atténuation dans la bande coupée de 40 dB. Spécifiez une fréquence d’échantillonnage de 2 000 Hz. Tracez la réponse en fréquence du filtre.

Fs = 2000;
[B,A] = ellip(20,0.1,40,[0.3 0.7],"ctf");
freqz(B,A,2048,Fs,"ctf")

Filtrez un signal chirp à balayage linéaire où la fréquence de Nyquist se produit à 1 seconde. Comparez les spectres des signaux d’entrée et de sortie.

t = 0:1/Fs:1;
x = chirp(t,0,t(end),Fs/2)';
y = filtfilt(B,A,x,"ctf");
pspectrum([x y],Fs,Leakage=1,FrequencyResolution=1)

Filtrage des artefacts sinusoïdaux dans un signal oscillatoire

Depuis R2024b

Recréez des artefacts sinusoïdaux bruités avec un filtrage à phase zéro par fonction de transfert. Filtrez un signal oscillatoire avec une CTF et des valeurs d’échelle.

Créez le signal u composé de bruit normalement distribué et de trois formes d’onde sinusoïdales. La fréquence d’échantillonnage est de 1 kHz.

rng("default")
Fs = 1e3;
ts = (0:1/Fs:2)';

a0 = [3 2 1];
f0 = [0.1 0.5 0.9]*Fs/2;
p0 = [0 pi/4 pi/2];

u = 0.1*randn(size(ts)) + 0.1*sin(2*pi*f0.*ts+p0)*a0';

Recréez des artefacts sinusoïdaux bruités en filtrant n0 avec un filtre de Butterworth numérique coupe-bande de troisième ordre et créez le signal v.

[b,a] = butter(3,[0.15 0.85],"stop");
v = filtfilt(b,a,u);

Comparez u et v. Remarquez que les deux signaux sont en phase.

tiledlayout("flow")
nexttile
strips([u(ts<0.1) v(ts<0.1)],0.1,Fs)
legend(["u" "v"],Location="southeast")
xlabel("Time (seconds)")
nexttile
pspectrum([u v],Fs)
legend(["u" "v"],Location="southeast")

Créez le signal oscillatoire commandé en tension x. Ajoutez les artefacts sinusoïdaux bruités représentés par le signal v.

vo = exp(-2*abs(ts-1)).*sin(8*pi*ts);
x = vco(vo,[0.25 0.75]*Fs/2,Fs) + v;

Filtrez le signal x avec un filtre de Tchebychev de type II d’ordre 24. Utilisez le format CTF et des valeurs d’échelle (B,A,g).

[B,A,g] = cheby2(24,50,[0.2 0.8],"ctf");
y = filtfilt({B,A,g},x);

Comparez l’amplitude au carré des transformées de Fourier à court terme. Remarquez la forte diminution d’amplitude dans les bandes atténuées.

tiledlayout("flow")
nexttile
stft(x,Fs,Window=bohmanwin(128),OverlapLength=120, ...
    FFTLength=512,FrequencyRange="onesided")
title("Input x")
nexttile
stft(y,Fs,Window=bohmanwin(128),OverlapLength=120, ...
    FFTLength=512,FrequencyRange="onesided")
title("Output y")

Arguments d'entrée

b, a — Coefficients de la fonction de transfert

vecteurs

Coefficients de la fonction de transfert, définis par des vecteurs. Si vous utilisez un filtre tout-pôle, définissez b à 1. Si vous utilisez un filtre tout-zéro, définissez a à 1.

Exemple : b = [1 3 3 1]/6 et a = [3 0 1 0]/3 définissent un filtre de Butterworth de troisième ordre avec une fréquence à 3 dB normalisée de 0,5π rad/échantillon.

Types de données : double | single

x — Signal d’entrée

vecteur | matrice | tableau ND

Signal d’entrée, défini par un vecteur, une matrice ou un tableau ND à valeurs réelles ou complexes. x doit être à valeurs finies. La longueur de x doit être supérieure au triple de l’ordre du filtre défini par max(length(B)-1, length(A)-1). La fonction s’applique le long de la première dimension de tableau de x sauf si x est un vecteur ligne. Si x est un vecteur ligne, la fonction s’applique le long de la deuxième dimension.

Exemple : cos(pi/4*(0:159))+randn(1,160) est un signal monocanal défini par un vecteur ligne.

Exemple : cos(pi./[4;2]*(0:159))'+randn(160,2) est un signal bicanal.

Types de données : double | single
Support des nombres complexes : Oui

sos — Coefficients des sections de second ordre

matrice

Coefficients des sections de second ordre, définis par une matrice. sos est une matrice L x 6 où le nombre de sections L doit être supérieur ou égal à 2. Si le nombre de sections est inférieur à 2, la fonction traite l’entrée comme un vecteur de numérateurs. Chaque ligne de sos correspond aux coefficients d’un filtre de second ordre (biquadratique). La i-ième ligne de sos correspond à [bi(1) bi(2) bi(3) ai(1) ai(2) ai(3)].

Exemple : s = [2 4 2 6 0 2;3 3 0 6 0 0] définit un filtre de Butterworth de troisième ordre avec une fréquence à 3 dB normalisée de 0,5π rad/échantillon.

Types de données : double | single

g — Valeurs d’échelle

scalaire | vecteur

Depuis R2024b

Valeurs d’échelle, définies par un scalaire à valeur réelle ou un vecteur à valeurs réelles à L + 1 éléments, où L est le nombre de sections CTF. Les valeurs d’échelle représentent la répartition du gain du filtre entre les sections de la représentation de filtre en cascade.

La fonction filtfilt applique un gain aux sections de filtre en utilisant la fonction scaleFilterSections selon la manière dont vous définissez g :

• Scalaire : la fonction répartit uniformément le gain entre toutes les sections de filtre.

• Vecteur : la fonction applique les L premières valeurs de gain aux sections de filtre correspondantes. La dernière valeur de gain est répartie uniformément entre toutes les sections de filtre.

Types de données : double | single

d — Filtre numérique

objet digitalFilter

Filtre numérique, défini par un objet digitalFilter. Utilisez designfilt pour générer un filtre numérique selon des spécifications de réponse en fréquence.

Exemple : d = designfilt("lowpassiir",FilterOrder=3,HalfPowerFrequency=0.5) définit un filtre de Butterworth de troisième ordre avec une fréquence à 3 dB normalisée de 0,5π rad/échantillon.

B,A — Coefficients de la CTF (fonction de transfert en cascade)

scalaires | vecteurs | matrices

Depuis R2024b

Coefficients de la fonction de transfert en cascade (CTF), définis par des scalaires, des vecteurs ou des matrices. B et A indiquent respectivement les coefficients de numérateur et de dénominateur de la fonction de transfert en cascade.

B doit être de taille L x (m + 1) et A doit être de taille L x (n + 1), où :

• L représente le nombre de sections de filtre.

• m représente l’ordre des numérateurs du filtre.

• n représente l’ordre des dénominateurs du filtre.

Pour plus d’informations sur le format de la fonction de transfert en cascade et les matrices de coefficients, veuillez consulter Spécifier des filtres numériques au format CTF.

Remarque

Si un élément de A(:,1) n’est pas égal à 1, la fonction filtfilt normalise les coefficients de filtre selon A(:,1). Dans ce cas, A(:,1) ne doit pas contenir de valeur nulle.

Types de données : double | single
Support des nombres complexes : Oui

Arguments de sortie

y — Signal filtré

vecteur | matrice | tableau ND

Signal filtré, renvoyé sous la forme d’un vecteur, d’une matrice ou d’un tableau ND.

• La fonction filtfilt renvoie y avec la même taille que x.

• Si vous spécifiez un argument d’entrée de type single, la fonction filtfilt calcule l’opération de filtrage en simple précision et renvoie y avec le type single. Sinon, filtfilt renvoie y avec le type double.

En savoir plus

Fonctions de transfert en cascade

Le partitionnement d’un filtre numérique IIR en plusieurs sections en cascade améliore sa stabilité numérique et réduit sa sensibilité aux erreurs de quantification des coefficients. Le format en cascade d’une fonction de transfert H(z) avec les L fonctions de transfert H₁(z), H₂(z), …, H_L(z) est le suivant :

$$H(z)={\displaystyle \prod _{l=1}^L{H}_l(z)}=H_1(z)\times H_2(z)\times \cdots \times H_L(z).$$

Spécifier des filtres numériques au format CTF

Vous pouvez spécifier des filtres numériques au format CTF pour l’analyse, la visualisation et le filtrage du signal. Pour spécifier un filtre, indiquez ses coefficients B et A. Vous pouvez également inclure le gain de mise à l’échelle du filtre pour toutes les sections en définissant g sous forme de scalaire ou de vecteur.

Coefficients des filtres

Lorsque vous définissez les coefficients sous forme de matrices à L lignes,

$$B=\left[\begin{array}{cccc}{b}_{11}& b_{12}& \cdots & b_{1,m+1}\\ b_{21}& b_{22}& \cdots & b_{2,m+1}\\ ⋮& ⋮& \ddots & ⋮\\ b_{L1}& b_{L2}& \cdots & b_{L,m+1}\end{array}\right],A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1,n+1}\\ a_{21}& a_{22}& \cdots & a_{2,n+1}\\ ⋮& ⋮& \ddots & ⋮\\ a_{L1}& a_{L2}& \cdots & a_{L,n+1}\end{array}\right],$$

le filtre est supposé avoir été spécifié comme une séquence de L fonctions de transfert en cascade. En conséquence, la fonction de transfert complète du filtre est la suivante :

$$H\left(z\right)=\frac{b_{11}+b_{12}{z}^{-1}+\cdots +b_{1,m+1}{z}^{-m}}{a_{11}+a_{12}{z}^{-1}+\cdots +a_{1,n+1}{z}^{-n}}\times \frac{b_{21}+b_{22}{z}^{-1}+\cdots +b_{2,m+1}{z}^{-m}}{a_{21}+a_{22}{z}^{-1}+\cdots +a_{2,n+1}{z}^{-n}}\times \cdots \times \frac{b_{L1}+b_{L2}{z}^{-1}+\cdots +b_{L,m+1}{z}^{-m}}{a_{L1}+a_{L2}{z}^{-1}+\cdots +a_{L,n+1}{z}^{-n}},$$

où est l’ordre du numérateur du filtre et n ≥ 0 est l’ordre du dénominateur.

• Si B et A sont tous les deux définis par des vecteurs, le système sous-jacent est supposé être un filtre IIR à section unique (L = 1), où B représente le numérateur de la fonction de transfert et A son dénominateur.

• Si B est un scalaire, le filtre est supposé être une cascade de filtres IIR tout-pôle dont chaque section a un gain total du système égal à B.

• Si A est un scalaire, le filtre est supposé être une cascade de filtres FIR dont chaque section a un gain total du système égal à 1/A.

Remarque

• Pour convertir des matrices de sections de second ordre en fonctions de transfert en cascade, utilisez la fonction sos2ctf.

• Pour convertir une représentation de filtre zéro-pôle-gain en fonctions de transfert en cascade, utilisez la fonction zp2ctf.

Coefficients et gain

Si vous utilisez un gain de mise à l’échelle total ou plusieurs gains de mise à l’échelle extraits des valeurs de coefficients, vous pouvez définir les coefficients et le gain avec un cell array au format {B,A,g}. La mise à l’échelle des sections de filtre est particulièrement importante lorsque vous utilisez l’arithmétique en virgule fixe, car elle permet de s’assurer que la sortie de chaque section a un niveau d’amplitude similaire. Vous évitez ainsi que la réponse du filtre comporte des inexactitudes dues à la précision numérique limitée.

Le gain peut être un scalaire représentant le gain total ou un vecteur de gains de section.

• Si le gain est un scalaire, sa valeur s’applique uniformément à toutes les sections du filtre en cascade.

• Si le gain est un vecteur, celui-ci doit contenir un élément de plus que le nombre de sections de filtre L dans la cascade. Chacune des L premières valeurs d’échelle s’applique à la section de filtre correspondante. La dernière valeur s’applique uniformément à toutes les sections du filtre en cascade.

Si vous définissez les matrices de coefficients et le vecteur de gains comme suit :

$$B=\left[\begin{array}{cccc}{b}_{11}& b_{12}& \cdots & b_{1,m+1}\\ b_{21}& b_{22}& \cdots & b_{2,m+1}\\ ⋮& ⋮& \ddots & ⋮\\ b_{L1}& b_{L2}& \cdots & b_{L,m+1}\end{array}\right],A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1,n+1}\\ a_{21}& a_{22}& \cdots & a_{2,n+1}\\ ⋮& ⋮& \ddots & ⋮\\ a_{L1}& a_{L2}& \cdots & a_{L,n+1}\end{array}\right],g=\left[\begin{array}{ccccc}{g}_1& g_2& \cdots & g_L& g_{\text{S}}\end{array}\right],$$

la fonction de transfert du système de filtre est supposée être la suivante :

$$H\left(z\right)=g_{\text{S}}\left(g_1\frac{b_{11}+b_{12}{z}^{-1}+\cdots +b_{1,m+1}{z}^{-m}}{a_{11}+a_{12}{z}^{-1}+\cdots +a_{1,n+1}{z}^{-n}}\times g_2\frac{b_{21}+b_{22}{z}^{-1}+\cdots +b_{2,m+1}{z}^{-m}}{a_{21}+a_{22}{z}^{-1}+\cdots +a_{2,n+1}{z}^{-n}}\times \cdots \times g_L\frac{b_{L1}+b_{L2}{z}^{-1}+\cdots +b_{L,m+1}{z}^{-m}}{a_{L1}+a_{L2}{z}^{-1}+\cdots +a_{L,n+1}{z}^{-n}}\right).$$

Conseils

• Vous pouvez obtenir des filtres au format CTF, y compris le gain de mise à l’échelle. Utilisez les sorties des fonctions de design de filtres IIR numériques telles que butter, cheby1, cheby2 et ellip. Indiquez l’argument de type de filtre "ctf" dans ces fonctions et spécifiez de renvoyer B, A et g pour obtenir les valeurs d’échelle. (depuis R2024b)

Références

[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] Lyons, Richard G. Understanding Digital Signal Processing. Upper Saddle River, NJ: Prentice Hall, 2004.

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

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

Capacités étendues

Génération de code C/C++

Générez du code C et C++ avec MATLAB® Coder™.

Notes d’usage et limitations :

Les objets digitalFilter ne sont pas supportés pour la génération de code.

Génération de code GPU

Générez du code CUDA® pour les GPU NVIDIA® avec GPU Coder™.

Notes d’usage et limitations :

Les objets digitalFilter ne sont pas supportés pour la génération de code.

GPU Arrays

Accélérez le code en exécutant les calculs sur une unité de traitement graphique (GPU) avec Parallel Computing Toolbox™.

Cette fonction supporte complètement les GPU arrays. Pour plus d’informations, veuillez consulter Exécuter les fonctions MATLAB sur un GPU (Parallel Computing Toolbox).

Historique des versions

Introduit avant R2006a

R2024b: Filtrer des signaux avec des fonctions de transfert en cascade

La fonction filtfilt supporte les entrées au format CTF (fonction de transfert en cascade).

Voir aussi

ctffilt | designfilt | digitalFilter | fftfilt | filter | filter2

Rubriques

• Remove the 60 Hz Hum from a Signal
• Practical Introduction to Digital Filtering
• Anti-Causal, Zero-Phase Filter Implementation