Est-il numériquement plus stable de mettre en œuvre le filtrage comme multiplication ou convolution?

J'écris un programme pour filtrer un signal de 20 000 échantillons avec un filtre Butterworth de cinquième ordre hors ligne. J'ai accès à une implémentation FFT. Il semble y avoir deux alternatives pour implémenter le filtrage:

• convolution du signal avec la réponse impulsionnelle dans le domaine temporel, ou
• multiplier le signal avec la réponse impulsionnelle dans le domaine fréquentiel et transformer le résultat en inverse

Ces méthodes seraient identiques dans le cas théorique du FT. Le faire dans la vraie vie avec la DFT, cependant, je suppose que les choses sont différentes. L'une des méthodes est-elle numériquement plus stable? Y a-t-il d'autres problèmes dont je devrais être conscient? Le nombre de calculs n'est pas important.

Avec la convolution, vous ne rencontrerez aucun problème de stabilité, car il n'y a pas de filtrage récursif, vous n'accumulerez donc aucune erreur. En d'autres termes, le système est entièrement composé de zéros, pas de pôles. J'ai entendu de façon anecdotique, mais non vérifiée par moi-même, que la convolution basée sur la FFT a une erreur plus faible que la convolution dans le domaine temporel, simplement parce qu'elle a des opérations arithmétiques O (n log n) plutôt que O (n ^ 2).

En règle générale, pour autant que je sache, les filtres Butterworth sont implémentés en tant que filtres récursifs (IIR), c'est donc un sujet différent. Les filtres IIR ont des pôles ainsi que des zéros, il peut donc y avoir des problèmes de stabilité dans la pratique. De plus, pour les filtres IIR, les méthodes basées sur la FFT ne sont pas une option, mais à la hausse, les filtres IIR ont tendance à être de très faible ordre.

En ce qui concerne les problèmes de stabilité avec les filtres IIR, ils ont tendance à avoir des problèmes à des ordres plus élevés - je vais simplement jeter un chiffre et dire que le 6e ordre le pousse. Au lieu de cela, ils sont généralement implémentés sous forme de biquads en cascade (sections de filtre de second ordre). Pour votre filtre du 5e ordre, écrivez-le comme une fonction de transfert de domaine z (ce sera une fonction rationnelle du 5e degré), puis factorisez-le dans ses 5 pôles et 5 zéros. Collectez les conjugués complexes et vous aurez deux biquades et un filtre de premier ordre. En général, les problèmes de stabilité ont tendance à apparaître lorsque les pôles se rapprochent du cercle unitaire.

Il peut également y avoir des problèmes avec les cycles de bruit et de limite dans les filtres IIR, il existe donc différentes topologies de filtre (c'est-à-dire forme directe I, forme directe II) qui ont des propriétés numériques différentes, mais je ne reviendrais pas sur ce point - utilisez simplement le double précision et ce sera certainement assez bon.

Dans presque tous les cas, votre meilleur choix n'est ni la convolution ni la FFT mais l'application directe du filtre IIR (en utilisant par exemple la fonction sosfilt ()). Ce sera beaucoup plus efficace en termes de consommation de CPU et de mémoire.

La différence numérique dépend du filtre spécifique. Le seul cas où une différence peut s'introduire est si les pôles sont très, très proches du cercle unitaire. Même là, quelques astuces peuvent vous aider. N'UTILISEZ PAS la représentation de la fonction de transfert et le filtre () mais utilisez des pôles et des zéros avec sosfilt (). Voici un exemple de différence.

n = 2^16;  % filter length
fs = 44100; % sample rate
x = zeros(n,1); x(1) = 1;
f0 = 15; % cutoff frequency in Hz
% design with poles and zeroes
[z,p,k] = butter(5,f0*2/fs);
clf
plot(sosfilt(zp2sos(z,p,k),x));
% design with transfer function
[b,a] = butter(5,f0*2/fs);
hold on
plot(filter(b,a,x),'k');

Le filtre () se détériore à une coupure d'environ 15 Hz à 44,1 kHz. Pour sosfilt (), la coupure peut être bien inférieure à 1/100 Hz à 44,1 kHz sans aucun problème.

SI vous avez des problèmes de stabilité, la FFT n'aide pas non plus. Puisque votre filtre est un filtre IIR, la réponse impulsionnelle est infinie et devrait être tronquée en premier. À ces fréquences très basses, la réponse impulsionnelle devient si longue que la FFT devient également impraticable.

Par exemple, si vous voulez une coupure de 1/100 Hz @ 44,1 kHz et que vous voulez une plage dynamique dans la réponse impulsionnelle de 100 dB, vous avez besoin d'environ 25 millions d'échantillons !!! C'est presque 10 minutes à 44,1 kHz et beaucoup, beaucoup plus longtemps que votre signal d'origine

Pourquoi pensez-vous que les choses seront différentes? Les concepts théoriques devraient se traduire par des applications pratiques, la seule différence étant celle des problèmes de virgule flottante, auxquels nous ne pouvons pas échapper. Vous pouvez facilement vérifier avec un exemple simple dans MATLAB:

x=randn(5,1);
y=randn(5,1);
X=fft(x,length(x)+length(y)-1);
Y=fft(y,length(x)+length(y)-1);

z1=conv(x,y);z2=ifft(X.*Y);
z1-z2

ans =

   1.0e-15 *

   -0.4441
   -0.6661
         0
   -0.2220
    0.8882
   -0.2220
         0
   -0.4441
    0.8882

Comme vous pouvez le constater, les erreurs sont de l'ordre de la précision de la machine. Il ne devrait y avoir aucune raison de ne pas utiliser la méthode FFT. Comme Endolith l'a mentionné, il est beaucoup plus courant d'utiliser l'approche FFT pour filtrer / convoluer / intercorréler, etc. et beaucoup plus rapidement, sauf pour de très petits échantillons (comme dans cet exemple). Non pas que le traitement dans le domaine temporel ne soit jamais fait ... tout se résume à l'application, aux besoins et aux contraintes.