Analog Discovery : étude d'un filtre

1. Introduction

Ce document montre comment utiliser un Analog Discovery (AD), avec l'interface Python développée dans Interface Python pour Analog Discovery, pour étudier un filtre linéaire.

Les scripts fonctionnent avec le fichier analog_0_1.py et le fichier dwfconstants.py (livré avec le SDK de l'AD).

2. Branchement du filtre

L'analog Discovery (version AD2 ou AD3) comporte deux sorties analogiques W1 et W2. Nous utilisons la sortie W1 pour fournir le signal appliqué à l'entrée du filtre. Il comporte aussi deux entrées analogiques différentiels 1+,1- et 2+,2-. Nous relions les bornes 1- et 2- à la masse et la borne 1+ à l'entrée du filtre, la borne 2+ à sa sortie. L'exemple traité est un circuit RLC dont voici le schéma :

Figure pleine page

La bobine est réalisée par un enroulement de 10 spires sur un noyau torique en ferrite. Son auto-inductance et sa résistance série à 10 kHz sont L=532 μF et r_b=0,37 Ω (valeurs données par un RLC-mètre). Le condensateur a une capacité C=979 nF et une résistance série r_c=0,15 Ω à 10 kHz. La résistance est R=10,0 Ω.

3. Test de linéarité

Le test de linéarité consiste à appliquer une tension sinusoïdale à l'entrée du filtre et à vérifier que la sortie est aussi sinusoïdale. En raison de l'impédance de la sortie analogique de l'AD, la tension V_e peut être différente de celle programmée.

analyseFiltreSinusoidal.py

import numpy as np
import matplotlib.pyplot as plt
from analog_0_1 import Device,AnalogInput,AnalogOutput
from scipy.signal.windows import get_window
from numpy.fft import fft
time sleep

La fonction calculSpectre calcule le spectre d'un signal par transformée de Fourier discrète (avec fenêtrage et ajout de zéros) :

def calculSpectre(s,Te,win,nz):
    # s : signal
    # Te : période d'échantillonnage 
    # win : fenêtrage
    # nz : nombre de blocs de zéros ajoutés
    N = len(s)
    if win=='rect':
        w = np.ones(N)
        a = 1.0
    elif win=='triang':
        w = get_window('triang',N)
        a = 0.5
    elif win=='hamming': 
        w = get_window('hamming',N)
        a = 0.54
    elif win=='hann':
        w = get_window('hann',N)
        a = 0.5
    elif win=='blackman':
        w = get_window('blackman',N)
        a = 0.42
    else:
        return None
    p = int(np.log(nz*N)/np.log(2))+1
    Np = 2**p
    u = np.zeros(Np)
    i1 = Np//2-(N-1)//2
    u[i1:i1+N] = s*w
    tfd = fft(u)*2/N/a
    freq = np.arange(Np)*1/(Np*Te)
    return freq,tfd

Une tension sinusoïdale de fréquence freq et d'amplitude de crête amp est programmée sur la sortie W1 puis on fait une acquisition sur les entrées 1 et 2 avec une fréquence d'échantillonnage 100 fois plus grande que la fréquence et une durée égale à 100 fois la période (soit un total de 10000 échantillons). Les tensions V_e(t) et V_s(t) sont tracées (sur une fenêtre de largeur égale à trois fois la période) ainsi que leur spectre en décibel.

freq = 5000
amp = 0.5
output.function(0,'sine',freq,amp)
output.start(0)
time.sleep(0.1)
fEchant = freq*100
duree = 100/freq
longueur = int(duree*fEchant)
time = analog.sampling(fEchant,longueur)
voltage = analog.record()
Ve = voltage[0,:]
Vs = voltage[1,:]
plt.figure()
plt.plot(time,Ve,label="Ve")
plt.plot(time,Vs,label="Vs")
plt.xlabel("t (s)",fontsize=16)
plt.ylabel("u (V)",fontsize=16)
plt.legend(loc="upper right",fontsize=16)
plt.xlim(0,4/freq)
plt.ylim(-1,1)
plt.grid()
Ve_eff = Ve.std()
Vs_eff = Vs.std()
valeurs = "Ve = %0.2f V, Vs = %0.2f V"%(Ve_eff,Vs_eff)
plt.title(valeurs)
print("Ve efficace = %f V"%Ve_eff)
print("Vs efficace = %f V"%Vs_eff)
Te = time[1]-time[0]
f,tfd_Ve = calculSpectre(Ve,Te,"blackman",8)
f,tfd_Vs = calculSpectre(Vs,Te,"blackman",8)
spectre_Ve = np.absolute(tfd_Ve)
spectre_Ve /= spectre_Ve.max()
spectre_Vs = np.absolute(tfd_Vs)
spectre_Vs /= spectre_Vs.max()
plt.figure()
plt.subplot(211)
plt.plot(f,20*np.log10(spectre_Ve),label="Ve")
plt.xlim(0,10*freq)
plt.ylabel("Ve",fontsize=16)
plt.grid()
plt.subplot(212)
plt.plot(f,20*np.log10(np.absolute(spectre_Vs)),label="Vs")
plt.xlabel("f (Hz)",fontsize=16)
plt.ylabel("Vs",fontsize=16)
plt.xlim(0,10*freq)
plt.grid()
plt.show()

Voici les résultats pour une fréquence de 7000 Hz (proche de la résonance) et une amplitude 0,5 V :

Le spectre de V_e révèle des harmoniques de rangs 2 et 3 qui sont absentes lorsque le filtre est débranché. Il s'agit donc d'un effet non linéaire (dû au noyau en ferrite). Cet effet peut être considéré comme négligeable puisque l'amplitude de l'harmonique le plus intense est à environ 50 dB du fondamental. Voici les courbes pour une amplitude de 0,8 V :

Pour cette amplitude, les non-linéarités ne sont plus négligeables. L'obtention de la réponse fréquentielle en régime linéaire devra donc se faire à un amplitude inférieure à 0,5 V.

4. Réponse fréquentielle

Il s'agit d'obtenir le gain et le déphasage en fonction de la fréquence pour une tension d'entrée sinusoïdale, dont l'amplitude est assez faible pour que le filtre fonctionne linéairement.

Posons

\begin{matrix} V_{e} (t) = E cos (ω t) & (1) \\ V_{s} (t) = S cos (ω t + φ) & (2) \end{matrix}

Si l'on note e_n et s_n les signaux échantillonnés correspondants, les valeurs moyennes sont définies par :

\begin{matrix} \bar{e} = \frac{1}{N} \sum_{n = 0}^{N - 1} e_{n} & (3) \\ \bar{s} = \frac{1}{N} \sum_{n = 0}^{N - 1} s_{n} & (4) \\ (5) \end{matrix}

La fonction numpy.mean permet de calculer la valeur moyenne.

L'amplitude efficace du signal retranchée de sa valeur moyenne est l'écart-type des N échantillons :

\begin{matrix} E_{e f f} = \sqrt{\frac{1}{N} \sum_{n = 0}^{N - 1} {(e_{n} - \bar{e})}^{2}} & (6) \\ S_{e f f} = \sqrt{\frac{1}{N} \sum_{n = 0}^{N - 1} {(s_{n} - \bar{s})}^{2}} & (7) \end{matrix}

La fonction numpy.std permet de calculer l'écart-type.

Le gain du filtre est le rapport des amplitudes efficaces :

G = \frac{S_{e f f}}{E_{e f f}} (8)

Afin de calculer le déphasage entre la sortie et l'entrée, on considère le produit suivant :

\begin{matrix} z (t) = & s (t) (V_{e} (t) - j V_{e} (t - \frac{T}{4})) & (9) \\ = & E S cos (ω t + φ) (cos (ω t) - j sin (ω t)) & (10) \\ = & \frac{E S}{2} (e^{j φ} + e^{- j (2 ω t + φ)}) & (11) \end{matrix}

Sa valeur moyenne est :

\bar{z} = \frac{E S}{2} e^{j φ} = E_{e f f} S_{e f f} e^{j φ} (12)

Le déphasage cherché est donc l'argument de $\bar{z}$ . On peut aussi calculer la valeur de la fonction de transfert :

\underset{̲}{H} = \frac{\bar{z}}{E_{e f f}^{2}} (13)

Le calcul des échantillons z_n nécessite de disposer du signal $e (t - \frac{T}{4})$ , c'est-à-dire le signal e(t) en avance d'un quart de période. Le nombre d'échantillons de décalage correspondant est la partie entière de N_p/4. Cependant, lorsque N_p est faible, l'obtention précise du signal en quadrature pose problème, à l'exception du cas peu probable où N_p est un nombre entier divisible par 4. Pour résoudre cette difficulté, nous effectuons une augmentation du nombre d'échantillons en procédant à une interpolation par transformée de Fourier. Cette méthode d'interpolation (d'un signal périodique) consiste à calculer la transformée de Fourier discrète du signal (TFD) à N echantillons, à séparer les deux parties conjuguées de taille N/2 puis à ajouter kN zéros entre ces deux parties. On obtient ainsi la TFD d'un signal qui comporte (k+1)N échantillons et qui est obtenu par la transformée de Fourier discrète inverse. Le nombre d'échantillons ajoutés par interpolation entre deux échantillons du signal initial est k.

analyseFiltreBode.py

import numpy as np
import matplotlib.pyplot as plt
from analog_0_1 import Device,AnalogInput,AnalogOutput
from scipy.signal.windows import get_window 
from numpy.fft import fft,ifft
from scipy.signal import firwin,convolve
import time

La fonction interpol effectue l'interpolation d'un signal échantillonné périodique par transformée de Fourier (ninter est le paramètre noté k ci-dessus) :

def interpol(x,ninter):
    # interpolation par FFT
    N = len(x)
    tfd = fft(x)
    N1 = N//2
    tfd2 = np.concatenate((tfd[0:N1],np.zeros(N*ninter),tfd[N1:N]))
    y = np.real(ifft(tfd2))*(ninter+1)
    return y

La fonction suivante calcule le déphasage :

def dephasage(t,u0,u1,ninter,freq):
    T = 1/freq
    te = t[1]-t[0]
    if ninter>0:
        v0 = interpol(u0,ninter)
        v1 = interpol(u1,ninter)
    d = int((ninter+1)*T/(4*te))
    N = len(v0)
    x = v1[d:N]*(v0[d:N]-1j*v0[0:N-d])
    phi = np.angle(x.mean())
    return phi

La fonction suivante effectue une mesure. Une tension sinusoïdale de fréquence freq et d'amplitude amp est programmée sur la sortie W1. Le paramètre duree désigne à la fois la durée de l'acquisition et la durée d'attente avant son démarrage. La fréquence d'échantillonnage est 100 fois la fréquence de la tension sinusoïdale appliquée.

def mesure(freq,amp,duree):
    output.function(0,'sine',freq,amp)
    output.start(0)
    time.sleep(duree)
    fEchant = freq*100
    longueur = int(duree*fEchant)
    t = analog.sampling(fEchant,longueur)
    voltage = analog.record()
    output.stop(0)
    Ve = voltage[0,:]
    Vs = voltage[1,:]
    Ve_eff = Ve.std()
    Vs_eff = Vs.std()
    G = Vs_eff/Ve_eff
    ninter = 4
    phi = dephasage(t,Ve,Vs,ninter,freq)
    return G,phi

Voici la boucle permettant de faire varier la fréquence :

f = np.logspace(2,5,100)
N = len(f)
list_G = np.zeros(N)
list_phi = np.zeros(N)
amp = 0.3
for k in range(N):
    duree = 80/f[k]
    G,phi = mesure(f[k],amp,duree)
    print(f[k],G,phi)
    list_G[k] = G
    list_phi[k] = phi
    
GdB = 20*np.log10(list_G)
np.savetxt("bode.txt",np.array([f,list_G,list_phi]).T)
plt.figure()
plt.subplot(211)
plt.plot(f,GdB,"b-")
plt.xlabel("f")
plt.ylabel("GdB")
plt.grid()
plt.xscale("log")
plt.subplot(212)
plt.plot(f,list_phi*180/np.pi,"b-")
plt.xlabel("f")
plt.ylabel("phi")
plt.grid()
plt.xscale("log")
plt.show()

La durée d'une acquisition est 80 fois la période de la sinusoïde et le nombre d'échantillons est donc 8000, un peu en dessous de la taille du tampon de l'Analog Discovery 2 (8192).

import numpy as np 
from matplotlib.pyplot import *
[f,G,phi] = np.loadtxt("bode-1.txt",unpack=True)
GdB = 20*np.log10(np.absolute(G))
figure(figsize=(12,8))
subplot(211)
plot(f,GdB,"b-")
ylabel("GdB",fontsize=16)
grid()
xscale("log")
subplot(212)
plot(f,phi*180/np.pi,"b-")
xlabel("f",fontsize=16)
ylabel("phi",fontsize=16)
grid()
xscale("log")

fig1.pdf

Afin de visualiser en détail le voisinage de la résonance, effectuons un balayage de fréquence linéaire sur un intervalle plus petit :

import numpy as np 
from matplotlib.pyplot import *
[f,G,phi] = np.loadtxt("bode-2.txt",unpack=True)
figure(figsize=(12,8))
subplot(211)
plot(f,G,"b-")
ylabel("G",fontsize=16)
ylim(0,1)
grid()
xscale("log")
subplot(212)
plot(f,phi*180/np.pi,"b-")
xlabel("f",fontsize=16)
ylabel("phi",fontsize=16)
grid()
xscale("log")

fig2.pdf

La fréquence de résonance est extraite :

i = np.argmax(G)
fr = f[i]
Gmax = G[i]

print(fr)
--> np.float64(6864.321608040202)

Nous pouvons tester le modèle présenté sur le schéma. Pour cela, il faut rechercher la valeur de r=r_b+r_c optimale. Étant donné le gain à la résonance (0,9), il semble pertinent de poser :

R = G_{m a x} (r + R)

La valeur précise de L est calculée à partir de la fréquence de résonance :

C = 979e-9
L = 1/((2*np.pi*fr)**2*C)
R = 10.0
r = (1-Gmax)/Gmax*R

print(L)
--> np.float64(0.0005491137586964649)

print(r)
--> np.float64(0.8801805007971174)

La valeur de L obtenue est notablement plus grande que celle donnée par le RLC-mètre.

f_modele = np.logspace(2,5,150)
w = f_modele*2*np.pi
H = R/(1j*L*w+r+1/(1j*C*w)+R)
G_modele = np.absolute(H)
GdB_modele = 20*np.log10(G_modele)
phi_modele = np.angle(H)
[f,G,phi] = np.loadtxt("bode-1.txt",unpack=True)
GdB = 20*np.log10(np.absolute(G))
figure(figsize=(12,8))
subplot(211)
plot(f,GdB,"b-")
plot(f_modele,GdB_modele,"k--")
ylabel("GdB",fontsize=16)
grid()
xscale("log")
subplot(212)
plot(f,phi*180/np.pi,"b-")
plot(f_modele,phi_modele*180/np.pi,"k--")
xlabel("f",fontsize=16)
ylabel("phi",fontsize=16)
grid()
xscale("log")

fig3.pdf

f_modele = np.linspace(4000,10000,200)
w = f_modele*2*np.pi
H = R/(1j*L*w+r+1/(1j*C*w)+R)
G_modele = np.absolute(H)
GdB_modele = 20*np.log10(G_modele)
phi_modele = np.angle(H)
[f,G,phi] = np.loadtxt("bode-2.txt",unpack=True)
figure(figsize=(12,8))
subplot(211)
plot(f,G,"b-")
plot(f_modele,G_modele,"k--")
ylabel("G",fontsize=16)
ylim(0,1)
grid()
xscale("log")
subplot(212)
plot(f,phi*180/np.pi,"b-")
plot(f_modele,phi_modele*180/np.pi,"k--")
xlabel("f",fontsize=16)
ylabel("phi",fontsize=16)
grid()
xscale("log")

fig4.pdf

L'accord entre le modèle est la réponse expérimentale est bon. On voit cependant que le gain du modèle est un peu inférieur à haute fréquence, ce qui peut être attribué à l'augmentation de r avec la fréquence. De même, le gain du modèle à basse fréquence est un peu plus grand, ce qui dû à une résistance r moindre. Un modèle plus précis devrait évidemment représenter les variations de r_b et de r_c avec la fréquence (la première augmente avec la fréquence, la seconde diminue). La valeur de r=r_b+r_c adoptée permet d'obtenir un très bon ajustement du modèle avec l'expérience au voisinage de la résonance. Il faut remarquer que l'expérience ne permet pas d'accéder à r_b et r_c séparément. D'après les mesures données par le RLC-mètre, ces deux résistances sont du même ordre de grandeur à la fréquence de résonance. Il semble que la variation de r avec la fréquence soit dominée par l'augmentation de r_b, qui est vraisemblablement plus grande que r_c.

5. Réponse indicielle

La réponse indicielle (ou réponse à un échelon) est obtenue en appliquant à l'entrée du filtre une tension nulle puis une tension positive.

On programme sur la sortie W1 un signal constitué d'une première moitié à zéro, la seconde moitié ayant la valeur amp. L'acquisition de l'entrée et de la sortie du filtre est faite avec un déclenchement sur le front montant de l'entrée. Par défaut le point de déclenchement est au milieu de l'enregistrement : la condition de déclenchement est donc effective seulement après une durée égale à la moitié de la durée de l'enregistrement. Pour réduire cette durée d'attente (et ainsi décaler le points de déclenchement vers la gauche), il faut attribuée une valeur positive à la durée définie par le paramètre position de la fonction analogTrigger.

analyseFiltreEchelon.py

import numpy as np
import matplotlib.pyplot as plt
from analog_0_1 import Device,AnalogInput,AnalogOutput

device = Device(-1)
device.open()
output = AnalogOutput(device)
analog = AnalogInput(device)
analog.channels([0,1],[5,5])

duree = 0.001
freq = 1/(2*duree)
Ne = 8000
fEchant = freq*Ne
samples = np.zeros(Ne)
samples[Ne//2:Ne] = 1
amp = 0.5
output.waveform(0,samples,freq,amp,periods=1,repeat=1)
longueur = int(2*duree*fEchant)
time = analog.sampling(fEchant,longueur)
analog.analogTrigger(0,amp/2,position=duree*0.5)
output.start(0)
voltage = analog.record()
output.stop(0)
Ve = voltage[0,:]
Vs = voltage[1,:]
plt.figure()
plt.plot(time*1e3,Ve,label="Ve")
plt.plot(time*1e3,Vs,label="Vs")
plt.xlabel("t (ms)",fontsize=16)
plt.ylabel("u (V)",fontsize=16)
plt.legend(loc="upper right",fontsize=16)
plt.ylim(-0.6,0.6)
plt.grid()
plt.show()

À la fin de la forme d'onde programmée, la sortie W1 est remise à zéro, ce qui permet d'observer aussi la réponse à un échelon négatif (marche descendante). Le même résultat pourrait être obtenu au moyen d'un signal périodique carré avec offset, ce que l'on fait généralement lorsqu'on utilise un générateur de signaux pour obtenir la réponse indicielle.