Intégrateur numérique

2.a. Conception du filtre

Soit y_n le signal numérique obtenu par échantillonnage de y(t). L'équation différentielle $\mathrm{(1)}$ peut aussi s'écrire :

$$\frac{dy\left(t\right)}{dt}=\frac{1}{\tau }x\left(t\right)\text{(3)}$$

En remplaçant la dérivée par une différence finie, on obtient :

$$y_n=y_{n-1}+\frac{T_e}{\tau }{x}_n\text{(4)}$$

La relation précédente est de la forme suivante :

$$a_0{y}_n=b_0{x}_n+b_1{x}_{n-1}-a_1{y}_{n-1}\text{(5)}$$

Cette relation de récurrence défini un filtre récursif (à réponse impulsionnelle infinie), dont la fonction de transfert en Z est :

$$H\left(z\right)=\frac{T_e}{\tau }\frac{1}{1-z^{-1}}\text{(6)}$$

Pour tracer sa réponse fréquentielle, on peut poser T_e/τ=1 puisque ce rapport n'a pas d'effet sur la forme du signal en sortie, seulement sur son amplitude :

import scipy.signal
import numpy
from matplotlib.pyplot import *

a=[1,-1]
b=[1]
w,h=scipy.signal.freqz(b,a)
figure()
subplot(211)
plot(w/(2*numpy.pi),20*numpy.log10(numpy.absolute(h)))
xlabel("f/fe")
ylabel("GdB")
xscale('log')
grid()
subplot(212)      
plot(w/(2*numpy.pi),numpy.angle(h)/numpy.pi)
xlabel("f/fe")
ylabel("phase/pi")
xscale('log')
grid()
              

figA.pdf

On n'obtient pas du tout le comportement intégrateur recherché, puisque le déphasage devrait être constamment égal à -π/2. La solution consiste à écrire un schéma numérique à deux pas (de type Adams-Moulton) :

$$y_n=y_{n-1}+\frac{T_e}{\tau }(x_n+x_{n-1})\text{(7)}$$

import scipy.signal
import numpy
from matplotlib.pyplot import *

a=[1,-1]
b=[1,1]
w,h=scipy.signal.freqz(b,a)
figure()
subplot(211)
plot(w/(2*numpy.pi),20*numpy.log10(numpy.absolute(h)))
xlabel("f/fe")
ylabel("GdB")
xscale('log')
grid()
subplot(212)      
plot(w/(2*numpy.pi),numpy.angle(h)/numpy.pi)
xlabel("f/fe")
ylabel("phase/pi")
xscale('log')
grid()
              

figB.pdf

On obtient ainsi un très bon intégrateur, sauf à proximité de la fréquence de Nyquist (f_e/2).

2.b. Réalisation du filtrage

Pour réaliser le filtrage d'une liste d'échantillons, on remarque que le calcul de la sortie commence à y₁, et qu'il faut choisir une valeur pour y₀. On prendra y₀=0.

Voici une fonction réalisant le filtrage d'une liste d'échantillons stockés en mémoire, pour un filtre récursif avec un numérateur et un dénominateur de degré 1 :

def filtrage_recursif(x,a,b):
    N = len(x)
    y = numpy.zeros(N)
    for n in range(1,N):
        y[n] = (-a[1]*y[n-1]+b[0]*x[n]+b[1]*x[n-1])/a[0]
    return y
                

On applique le filtre précédent à une impulsion :

x = numpy.zeros(100)
x[0] = 1.0
y = filtrage_recursif(x,a,b)
figure()
stem(y)
xlabel("n")
ylabel("y(n)")
axis([0,100,0,2])
grid()
                

figC.pdf

La réponse impulsionnelle d'un intégrateur est un échelon. Elle ne tend pas vers zéro lorsque n tend vers l'infini (elle est constante), ce qui montre que l'intégrateur est instable.

2.c. Test de l'intégrateur

Un bon moyen de tester un intégrateur est de lui fournir un signal créneau. Nous allons donc numériser un signal créneau délivré par un GBF avec la carte d'acquisition et le script suivant :

# -*- coding: utf-8 -*-
import pycan.main as pycan
import numpy
from matplotlib.pyplot import *
import scipy.signal
fe=40000.0
T=0.1
te=1.0/fe
N = int(T/te)
can = pycan.Sysam("SP5")
can.config_entrees([0],[10.0])
can.config_echantillon(te*10**6,N)
can.acquerir()
t0=can.temps()[0]
u0=can.entrees()[0]
can.fermer()
numpy.savetxt("signal-1.txt",[t0,u0])
            

On applique le filtre intégrateur. Le gain est réduit en choisissant une valeur plus faible de b₀=b₁.

[t0,u0] = numpy.loadtxt("signal-1.txt")
a=[1,-1]
b=[0.01,0.01]
y = filtrage_recursif(u0,a,b)
figure()
plot(t0,u0)
plot(t0,y)
xlabel("t (s)")
axis([0,0.05,-2,2])
grid()
            

figD.pdf

On observe une forte dérive, due à la présence d'une composante de fréquence nulle dans le signal x_n, bien que celle-ci soit très faible.

Lorsque le signal x(t) comporte une composante de fréquence nulle que l'on souhaite intégrer, il faut utiliser cet intégrateur. Il faut cependant faire une compensation qui permette d'annuler tout décalage accidentel.

3.a. Filtre analogique

Pour stabiliser le filtre intégrateur afin de ne pas avoir de dérive en sortie, il faut ramener le gain à fréquence nulle à une valeur finie. Un moyen simple de le faire est d'utiliser un filtre passe-bas du premier ordre.

Considérons un filtre passe-bas dont la fonction de transfert est :

$$H\left(\omega \right)=\frac{g}{1+j\tau \omega }\text{(8)}$$

Ce filtre est quasiment intégrateur lorsque $\omega \gg 1/\tau$.

L'équation différentielle associée à ce filtre est :

$$y\left(t\right)+\tau \frac{dy\left(t\right)}{dt}=gx\left(t\right)\text{(9)}$$

3.b. Filtre numérique

La discrétisation de l'équation différentielle conduit à la relation suivante :

$$y_n+\tau \frac{y_n-y_{n-1}}{T_e}=g\frac{x_n+x_{n-1}}{2}\text{(10)}$$

On obtient ainsi un filtre récursif défini par :

$$y_n=\frac{1}{1+\frac{T_e}{\tau }}{y}_{n-1}+\frac{g}{2}\frac{\frac{T_e}{\tau }}{1+\frac{T_e}{\tau }}(x_n+x_{n-1})\text{(11)}$$

Sa fonction de transfert en Z est :

$$H\left(z\right)=\frac{g}{2}r\frac{T_e}{\tau }\frac{1+z^{-1}}{1-r{z}^{-1}}\text{(12)}$$

avec :

$$r=\frac{1}{1+\frac{T_e}{\tau }}\text{(13)}$$

Le pôle est z=r qui est strictement inférieur à 1, ce qui rend le filtre stable. Pour que le domaine de fréquence d'intégration soit large, il faut que le rapport T_e/τ soit faible, et donc que r soit proche de 1. Par exemple, si l'on veut un domaine d'intégration qui commence à environ un centième de la fréquence d'échantillonnage, on choisit T_e/τ=1/100 ou moins. Pour la forme de la réponse fréquentielle, on peut choisir librement le coefficient en facteur dans H, qui n'a d'effet que sur l'amplitude du signal en sortie, pas sur sa forme.

rt = 0.005
r=1.0/(1+rt)
a=[1,-r]
b=[1,1]
w,h=scipy.signal.freqz(b,a,worN=numpy.logspace(-4,-0.31,1000)*2*numpy.pi)
figure()
subplot(211)
plot(w/(2*numpy.pi),20*numpy.log10(numpy.absolute(h)))
xlabel("f/fe")
ylabel("GdB")
xscale('log')
grid()
subplot(212)      
plot(w/(2*numpy.pi),numpy.angle(h)/numpy.pi)
xlabel("f/fe")
ylabel("phase/pi")
xscale('log')
grid()

                 

figE.pdf

Ce filtre passe-bas a une fréquence de coupure environ égale au millième de la fréquence d'échantillonnage, et un domaine intégrateur qui commence environ au centième de cette fréquence. Voyons la réponse impulsionnelle de ce filtre :

x = numpy.zeros(1000)
x[0] = 1.0
y = filtrage_recursif(x,a,b)
figure()
plot(y,".")
xlabel("n")
ylabel("y(n)")
axis([0,1000,0,2])
grid()
                

figF.pdf

La réponse impulsionnelle tend vers zéro car le filtre est stable, mais la décroissance est lente, d'autant plus que r est proche de 1. Le temps de réponse est de l'ordre de τ. La contrepartie d'un filtre intégrateur à très large bande est donc un temps de réponse très long, qui peut être gênant lors des transitoires (changement de fréquence par exemple). Il faudra donc choisir τ au plus juste en fonction de la fréquence minimale des signaux à intégrer.

3.c. Test du filtre

Voici l'application du filtre au signal créneau. Les valeurs de b₀=b₁ sont choisie pour abaisser le gain.

b=[0.01,0.01]
y = filtrage_recursif(u0,a,b)
figure()
plot(t0,u0)
plot(t0,y)
xlabel("t (s)")
axis([0,0.05,-2,2])
grid()
            

figG.pdf

Après le régime transitoire de durée τ, on observe bien un état stationnaire avec un décalage en sortie relativement important mais constant.

Pour obtenir une intégration vraie, correspondant à τ=1 dans la relation $\mathrm{(1)}$, il faut utiliser :

$$y_n=r{y}_{n-1}+\frac{T_e}{2}(x_n+x_{n-1})\text{(14)}$$

te = t0[1]-t0[0]
b=[te/2,te/2]
y = filtrage_recursif(u0,a,b)
figure()
plot(t0,y)
xlabel("t (s)")
axis([0,0.05,-1e-3,1e-3])
grid()
            

figH.pdf