2. Exemples🔗
2.1. Mean power (fonction, argparse)🔗
1#!/usr/bin/env python3
2
3"""
4Exemple de script (shebang, docstring, etc.) permettant une
5utilisation en module (`import mean_power`) et en exécutable (`python
6mean_power.py -h`);
7"""
8
9
10def mean_power(alist, power=1):
11 r"""
12 Retourne la racine `power` de la moyenne des éléments de `alist` à
13 la puissance `power`:
14
15 .. math:: \mu = (\frac{1}{N}\sum_{i=0}^{N-1} x_i^p)^{1/p}
16
17 `power=1` correspond à la moyenne arithmétique, `power=2` au *Root
18 Mean Squared*, etc.
19
20 Exemples:
21 >>> mean_power([1, 2, 3])
22 2.0
23 >>> mean_power([1, 2, 3], power=2)
24 2.160246899469287
25 """
26
27 # *mean* = (somme valeurs**power / nb valeurs)**(1/power)
28 mean = (sum( val ** power for val in alist ) / len(alist)) ** (1 / power)
29
30 return mean
31
32
33if __name__ == '__main__':
34
35 # start-argparse
36 import argparse
37
38 parser = argparse.ArgumentParser()
39 parser.add_argument('list', nargs='*', type=float, metavar='nombres',
40 help="Liste de nombres à moyenner")
41 parser.add_argument('-i', '--input',
42 nargs='?', type=argparse.FileType('r'),
43 help="Fichier contenant les nombres à moyenner")
44 parser.add_argument('-p', '--power',
45 type=float, default=1.,
46 help="'Puissance' de la moyenne (par défaut: %(default)s)")
47
48 args = parser.parse_args()
49 # end-argparse
50
51 if args.input: # Lecture des coordonnées du fichier d'entrée
52 # Le fichier a déjà été ouvert en lecture par argparse (type=file)
53 try:
54 args.list = [float(row) for row in args.input
55 if not row.strip().startswith('#')]
56 except ValueError:
57 parser.error("Impossible de déchiffrer le fichier {args.input!r}")
58
59 # Vérifie qu'il y a au moins un nombre dans la liste
60 if not args.list:
61 parser.error("La liste doit contenir au moins un nombre")
62
63 # Calcul
64 moyenne = mean_power(args.list, args.power)
65
66 # Affichage du résultat
67 print(f"La moyenne puissance 1/{args.power} des {len(args.list)} nombres "
68 f"à la puissance {args.power} est {moyenne}.")
Source: mean_power.py
2.2. Animal (POO)🔗
1#!/usr/bin/env python3
2
3"""
4Exemple (tragique) de Programmation Orientée Objet.
5"""
6
7
8# Définition d'une classe ==============================
9
10class Animal:
11 """
12 Un animal, défini par sa `masse` et sa vitalité (`estVivant`).
13 """
14
15 def __init__(self, masse):
16 """
17 Initialisation d'un Animal, a priori vivant.
18
19 :param float masse: masse en kg (> 0)
20 :raise ValueError: masse non réelle ou négative
21 """
22
23 self.masse = float(masse) # Ajout d'un attribut "masse"
24 if self.masse < 0:
25 raise ValueError("La masse ne peut pas être négative.")
26
27 self.estVivant = True # Ajout d'un attribut "estVivant"
28
29
30 def __str__(self):
31 """
32 Surcharge de la fonction `str()`.
33
34 L'affichage *informel* de l'objet dans l'interpréteur.
35
36 P.ex. `print(a)` sera résolu comme `a.__str__()`, ou `a.__repr__()` en
37 l'absence de `__str__`.
38
39 :return: une chaîne de caractères
40 """
41
42 return f"Animal {'vivant' if self.estVivant else 'mort'}, " \
43 f"{self.masse:.0f} kg"
44
45
46 def meurt(self):
47 """
48 L'animal meurt.
49 """
50
51 self.estVivant = False
52
53
54 def grossit(self, masse):
55 """
56 L'animal grossit (ou maigrit) d'une certaine masse (valeur algébrique).
57
58 :param float masse: prise (>0) ou perte (<0) de masse.
59 :raise ValueError: masse non réelle.
60 """
61
62 self.masse += float(masse)
63
64
65 def __repr__(self):
66 """
67 *AVANCÉ:* Surcharge de la fonction `repr()`.
68
69 L'affichage *formel* de l'objet dans l'interpréteur.
70
71 :return: une chaîne de caractères
72 """
73
74 return f"{self.__class__.__name__}({self.masse})"
75
76
77# Définition d'une classe héritée ==============================
78
79class AnimalFeroce(Animal):
80 """
81 Un animal féroce est un animal qui peut dévorer d'autres animaux.
82
83 La classe-fille hérite des attributs et méthodes de la
84 classe-mère, mais peut les surcharger (i.e. en changer la
85 définition), ou en ajouter de nouveaux:
86
87 - la méthode `AnimalFeroce.__init__()` dérive directement de
88 `Animal.__init__()` (même méthode d'initialisation);
89 - `AnimalFeroce.__str__()` surcharge `Animal.__str__()`;
90 - `AnimalFeroce.devorer()` est une nouvelle méthode propre à
91 `AnimalFeroce`.
92 """
93
94 def __str__(self):
95 """
96 Surcharge de la fonction `str()`.
97 """
98
99 return "Animal féroce " \
100 f"{'bien vivant' if self.estVivant else 'mais mort'}, " \
101 f"{self.masse:.0f} kg"
102
103 def devore(self, other):
104 """
105 L'animal (self) devore un autre animal (other).
106
107 * Si other est également un animal féroce, il faut que self soit plus
108 gros que other pour le dévorer. Sinon, other se défend et self meurt.
109 * Si self dévore other, other meurt, self grossit de la masse de other
110 (jusqu'à 10% de sa propre masse) et other maigrit d'autant.
111
112 :param Animal other: animal à dévorer
113 :return: prise de masse (0 si self meurt)
114 """
115
116 if isinstance(other, AnimalFeroce) and (other.masse > self.masse):
117 # Pas de chance...
118 self.meurt()
119 prise = 0.
120 else:
121 other.meurt() # Other meurt
122 prise = min(other.masse, self.masse * 0.1)
123 self.grossit(prise) # Self grossit
124 other.grossit(-prise) # Other maigrit
125
126 return prise
127
128
129# Définition d'une autre classe héritée ==============================
130
131class AnimalGentil(Animal):
132 """
133 Un animal gentil est un animal avec un `nom`.
134
135 La classe-fille hérite des attributs et méthodes de la
136 classe-mère, mais peut les surcharger (i.e. en changer la
137 définition), ou en ajouter de nouveaux:
138
139 - la méthode `AnimalGentil.__init__()` surcharge l'initialisation originale
140 `Animal.__init__()`;
141 - `AnimalGentil.__str__()` surcharge `Animal.__str__()`;
142 """
143
144 def __init__(self, masse, nom='Youki'):
145 """
146 Initialisation d'un animal gentil, avec son masse et son nom.
147 """
148
149 # Initialisation de la classe parente (nécessaire pour assurer
150 # l'héritage)
151 Animal.__init__(self, masse)
152
153 # Attributs propres à la classe AnimalGentil
154 self.nom = nom
155
156 def __str__(self):
157 """
158 Surcharge de la fonction `str()`.
159 """
160
161 return f"{self.nom}, un animal gentil " \
162 f"{'bien vivant' if self.estVivant else 'mais mort'}, " \
163 f"{self.masse:.0f} kg"
164
165 def meurt(self):
166 """
167 L'animal gentil meurt, avec un éloge funéraire.
168 """
169
170 Animal.meurt(self)
171 print(f"Pauvre {self.nom} meurt, paix à son âme...")
172
173
174if __name__ == '__main__':
175
176 # Exemple d'utilisation des classes définies ci-dessus
177
178 print("Une tragédie en trois actes".center(70, '='))
179
180 print("Acte I: la vache prend 10 kg.".center(70, '-'))
181 vache = Animal(500.) # Instantiation d'un animal de 500 kg
182 vache.grossit(10) # La vache grossit de 10 kg
183 print(vache)
184
185 print("Acte II: Dumbo l'éléphant".center(70, '-'))
186 elephant = AnimalGentil(1000., "Dumbo") # Instantiation d'un animal gentil
187 print(elephant)
188
189 print("Acte III: le féroce lion".center(70, '-'))
190 lion = AnimalFeroce(200) # Instantiation d'un animal féroce
191 print(lion)
192
193 print("Scène tragique: le lion dévore l'éléphant...".center(70, '-'))
194 lion.devore(elephant) # Le lion dévore l'éléphant
195
196 print(elephant)
197 print(lion)
Exécution du code:
$ ./animal.py
=====================Une tragédie en trois actes======================
--------------------Acte I: la vache prend 10 kg.---------------------
Animal vivant, 510 kg
----------------------Acte II: Dumbo l'éléphant-----------------------
Dumbo, un animal gentil bien vivant, 1000 kg
-----------------------Acte III: le féroce lion-----------------------
Animal féroce bien vivant, 200 kg
-------------Scène tragique: le lion dévore l'éléphant...-------------
Pauvre Dumbo, paix à son âme...
Dumbo, un animal gentil mais mort, 980 kg
Animal féroce bien vivant, 220 kg
Source: animal.py
2.3. Cercle circonscrit (POO, argparse)🔗
1#!/usr/bin/env python3
2
3"""
4Calcule le cercle circonscrit à 3 points du plan.
5
6Ce script sert d'illustration à plusieurs concepts indépendants:
7
8- un exemple de script (shebang, docstring, etc.) permettant une
9 utilisation en module (`import circonscrit`) et en exécutable
10 (`python circonscrit.py -h`);
11- des exemples de Programmation Orientée Objet: classe `Point` et la
12 classe héritière `Vector`;
13- un exemple d'utilisation du module `argparse` de la bibliothèque
14 standard, permettant la gestion des arguments de la ligne de
15 commande;
16- l'utilisation de tests unitaires sous la forme de `doctest` (tests
17 inclus dans les *docstrings* des éléments à tester).
18
19 Pour exécuter les tests unitaires du module:
20
21 - avec doctest: `python -m doctest -v circonscrit.py`;
22 - avec pytests: `py.test --doctest-modules -v circonscrit.py`;
23"""
24
25__author__ = "Yannick Copin <y.copin@ipnl.in2p3.fr>"
26__version__ = "Time-stamp: <2014-01-12 22:19 ycopin@lyonovae03.in2p3.fr>"
27
28# Définition d'une classe ==============================
29
30
31class Point:
32
33 """
34 Classe définissant un `Point` du plan, caractérisé par ses
35 coordonnées `x`,`y`.
36 """
37
38 def __init__(self, x, y):
39 """
40 Méthode d'instanciation à partir de deux coordonnées réelles.
41
42 >>> Point(0, 1) # doctest: +ELLIPSIS
43 <....Point object at 0x...>
44 >>> Point(1 + 3j)
45 Traceback (most recent call last):
46 ...
47 TypeError: Point.__init__() missing 1 required positional argument: 'y'
48 """
49
50 try: # Convertit les coords en `float`
51 self.x = float(x)
52 self.y = float(y)
53 except (ValueError, TypeError) as exc:
54 raise TypeError(f"Invalid input coordinates ({x}, {y})") from exc
55
56 def __str__(self):
57 """
58 Surcharge de la fonction `str()`: l'affichage *informel* de l'objet
59 dans l'interpréteur, p.ex. `str(self)` sera résolu comme
60 `self.__str__()`
61
62 Retourne une chaîne de caractères.
63
64 >>> str(Point(1, 2))
65 'Point (x=1.0, y=2.0)'
66 """
67
68 return f"Point (x={self.x}, y={self.y})"
69
70 def is_origin(self):
71 """
72 Teste si le point est à l'origine en testant la nullité des deux
73 coordonnées.
74
75 Attention aux éventuelles erreurs d'arrondis: il faut tester
76 la nullité à la précision numérique près.
77
78 >>> Point(0, 0).is_origin()
79 True
80 >>> Point(1, 2).is_origin()
81 False
82 """
83
84 from sys import float_info
85
86 # float_info.epsilon est le plus petit float non nul
87
88 return ((abs(self.x) <= float_info.epsilon) and
89 (abs(self.y) <= float_info.epsilon))
90
91 def distance(self, other):
92 """
93 Méthode de calcul de la distance du point (`self`) à un autre point
94 (`other`).
95
96 >>> A = Point(1,0); B = Point(1,1); A.distance(B)
97 1.0
98 """
99
100 # hypot(dx, dy) = sqrt(dx**2 + dy**2)
101 return ((self.x - other.x) ** 2 + (self.y - other.y) ** 2) ** 0.5
102
103
104# Définition du point origine O
105O = Point(0, 0)
106
107
108# Héritage de classe ==============================
109
110
111class Vector(Point):
112
113 """
114 Un `Vector` hérite de `Point` avec des méthodes additionnelles
115 (p.ex. la négation d'un vecteur, l'addition de deux vecteurs, ou
116 la rotation d'un vecteur).
117 """
118
119 def __init__(self, A, B):
120 """
121 Définit le vecteur `AB` à partir des deux points `A` et `B`.
122
123 >>> Vector(Point(1, 0), Point(1, 1)) # doctest: +ELLIPSIS
124 <....Vector object at 0x...>
125 >>> Vector(0, 1)
126 Traceback (most recent call last):
127 ...
128 AttributeError: 'int' object has no attribute 'x'
129 """
130
131 # Initialisation de la classe parente
132 Point.__init__(self, B.x - A.x, B.y - A.y)
133
134 # Attribut propre à la classe dérivée
135 self.sqnorm = self.x**2 + self.y**2 # Norme du vecteur au carré
136
137 def __str__(self):
138 """
139 Surcharge de la fonction `str()`: `str(self)` sera résolu comme
140 `Vector.__str__(self)` (et non pas comme `Point.__str__(self)`)
141
142 >>> str(Vector(O, Point(0, 1)))
143 'Vector (x=0.0, y=1.0)'
144 """
145
146 return f"Vector (x={self.x}, y={self.y})"
147
148 def __eq__(self, other):
149 """
150 Surcharge du test d'égalité `{self}=={other}`: l'instruction sera
151 résolue comme `self.__eq__(other)`.
152
153 >>> Vector(O, Point(0, 1)) == Vector(Point(1, 0), Point(1, 1))
154 True
155 """
156
157 # On teste ici la nullité de la différence des 2 vecteurs.
158 # D'autres tests auraient été possibles -- égalité des
159 # coordonnées, nullité de la norme de la différence, etc. --
160 # mais on tire profit de la méthode héritée `Point.is_origin()`
161 # testant la nullité des coordonnées (à la précision numérique
162 # près).
163 return (self - other).is_origin()
164
165 def __add__(self, other):
166 """
167 Surcharge de l'opérateur binaire `{self} + {other}`: l'instruction
168 sera résolue comme `self.__add__(other)`.
169
170 On construit une nouvelle instance de `Vector` à partir des
171 coordonnées propres à l'objet `self` et à l'autre opérande
172 `other`.
173
174 >>> A = Point(1, 0); B = Point(1, 1)
175 >>> Vector(A, B) + Vector(B, O) == Vector(A, O)
176 True
177 """
178
179 return Vector(O, Point(self.x + other.x, self.y + other.y))
180
181 def __sub__(self, other):
182 """
183 Surcharge de l'opérateur binaire `{self} - {other}`: l'instruction
184 sera résolue comme `self.__sub__(other)`.
185
186 Attention: ne surcharge pas l'opérateur unaire `-{self}`, géré
187 par `__neg__`.
188
189 >>> A = Point(1, 0); B = Point(1, 1)
190 >>> str(Vector(A, B) - Vector(A, B)) # Différence
191 'Vector (x=0.0, y=0.0)'
192 >>> -Vector(A, B) # Négation
193 Traceback (most recent call last):
194 ...
195 TypeError: bad operand type for unary -: 'Vector'
196 """
197
198 return Vector(O, Point(self.x - other.x, self.y - other.y))
199
200 def __abs__(self):
201 """
202 Surcharge la fonction `abs()` pour retourner la norme du vecteur.
203
204 >>> abs(Vector(O, Point(1, 0)))
205 1.0
206 """
207
208 # On pourrait utiliser sqrt(self.sqnorm), mais c'est pour
209 # illustrer l'utilisation de la méthode héritée
210 # `Point.distance`...
211 return Point.distance(self, O)
212
213 def rotate(self, angle, deg=False):
214 """
215 Rotation (dans le sens trigonométrique) du vecteur par un `angle`,
216 exprimé en radians ou en degrés.
217
218 >>> Vector(O, Point(0, 1)).rotate(90, deg=True) == Vector(O, Point(-1, 0))
219 True
220 """
221
222 from cmath import rect # Bibliothèque de fonctions complexes
223
224 # On calcule la rotation en passant dans le plan complexe
225 z = complex(self.x, self.y)
226 phase = angle if not deg else angle / 57.29577951308232 # [rad]
227 zu = z * rect(1.0, phase) # Rotation complexe (* exp(i*phase))
228
229 return Vector(O, Point(zu.real, zu.imag))
230
231
232def circumscribedCircle(M, N, P):
233 """
234 Calcule le centre et le rayon du cercle circonscrit aux points
235 M, N, P.
236
237 Retourne: (centre [Point], rayon [float])
238
239 Lève une exception `ValueError` si le rayon ou le centre du cercle
240 circonscrit n'est pas défini.
241
242 >>> M = Point(-1, 0); N = Point(1, 0); P = Point(0, 1)
243 >>> C, r = circumscribedCircle(M, N, P) # Centre O, rayon 1
244 >>> C.distance(O), round(r, 6)
245 (0.0, 1.0)
246 >>> circumscribedCircle(M, O, N) # Indéfini
247 Traceback (most recent call last):
248 ...
249 ValueError: Undefined circumscribed circle radius.
250 """
251
252 MN = Vector(M, N)
253 NP = Vector(N, P)
254 PM = Vector(P, M)
255
256 # Rayon du cercle circonscrit
257 m = abs(NP) # |NP|
258 n = abs(PM) # |PM|
259 p = abs(MN) # |MN|
260
261 d = (m + n + p) * (-m + n + p) * (m - n + p) * (m + n - p)
262 if d > 0:
263 rad = m * n * p / d**0.5
264 else:
265 raise ValueError("Undefined circumscribed circle radius.")
266
267 # Centre du cercle circonscrit
268 d = -2 * (M.x * NP.y + N.x * PM.y + P.x * MN.y)
269 if d == 0:
270 raise ValueError("Undefined circumscribed circle center.")
271
272 om2 = Vector(O, M).sqnorm # |OM|**2
273 on2 = Vector(O, N).sqnorm # |ON|**2
274 op2 = Vector(O, P).sqnorm # |OP|**2
275
276 x0 = -(om2 * NP.y + on2 * PM.y + op2 * MN.y) / d
277 y0 = (om2 * NP.x + on2 * PM.x + op2 * MN.x) / d
278
279 return (Point(x0, y0), rad) # (centre [Point], R [float])
280
281
282if __name__ == "__main__":
283
284 # start-argparse
285 import argparse
286
287 parser = argparse.ArgumentParser(
288 usage="%(prog)s [-p/--plot] [-i/--input coordfile | x1,y1 x2,y2 x3,y3]",
289 description="Compute the circumscribed circle to 3 points in the plan.",
290 )
291 parser.add_argument(
292 "coords",
293 nargs="*", type=str, metavar="x,y",
294 help="Coordinates of point"
295 )
296 parser.add_argument(
297 "-i", "--input",
298 nargs="?", type=argparse.FileType("r"),
299 help="Coordinate file (one 'x,y' per line)",
300 )
301 parser.add_argument(
302 "-P", "--plot",
303 action="store_true", default=False,
304 help="Draw the circumscribed circle",
305 )
306 parser.add_argument(
307 "-T", "--tests",
308 action="store_true", default=False,
309 help="Run doc tests"
310 )
311 parser.add_argument("--version", action="version", version=__version__)
312
313 args = parser.parse_args()
314 # end-argparse
315
316 if args.tests: # Auto-test mode
317 import sys
318 import doctest
319
320 fails, tests = doctest.testmod(verbose=True) # Run doc tests
321 sys.exit(fails > 0)
322
323 if args.input: # Lecture des coordonnées du fichier d'entrée
324 # Le fichier a déjà été ouvert en lecture par argparse (type=file)
325 args.coords = [coords
326 for coords in args.input
327 if not coords.strip().startswith("#")]
328
329 if len(args.coords) != 3: # Vérifie le nb de points
330 parser.error("Specify 3 points by their coordinates 'x,y' "
331 f"(got {len(args.coords)})")
332
333 points = [] # Liste des points
334 for i, arg in enumerate(args.coords, start=1):
335 try: # Déchiffrage de l'argument 'x,y'
336 x, y = (float(t) for t in arg.split(","))
337 except ValueError:
338 parser.error(f"Cannot decipher coordinates #{i}: {arg!r}")
339
340 points.append(Point(x, y)) # Création du point et ajout à la liste
341 print(f"#{i:d}: {points[-1]}") # Affichage du dernier point
342
343 # Calcul du cercle cisconscrit (lève une ValueError en cas de problème)
344 center, radius = circumscribedCircle(*points) # Délistage
345 print(f"Circumscribed circle: {center}, radius: {radius}")
346
347 if args.plot: # Figure
348 import matplotlib.pyplot as plt
349
350 fig, ax = plt.subplots()
351 ax.set(aspect="equal")
352 # Points
353 ax.plot([p.x for p in points], [p.y for p in points], "ko")
354 for i, p in enumerate(points, start=1):
355 ax.annotate(f"#{i}", (p.x, p.y),
356 xytext=(5, 5), textcoords="offset points")
357 # Cercle circonscrit
358 c = plt.matplotlib.patches.Circle((center.x, center.y),
359 radius=radius, fc="none", ec="k")
360 ax.add_patch(c) # Cercle
361 ax.plot(center.x, center.y, "r+") # Centre
362
363 plt.show()
$ ./circonscrit.py -h
usage: circonscrit.py [-p/--plot] [-i/--input coordfile | x1,y1 x2,y2 x3,y3]
positional arguments:
x,y Coordinates of point
optional arguments:
-h, --help show this help message and exit
-i [INPUT], --input [INPUT]
Coordinate file (one 'x,y' per line)
-p, --plot Draw the circumscribed circle
--version show program's version number and exit
$ ./circonscrit.py -p 0,0 1,0 1,1
Source: circonscrit.py
2.4. Matplotlib🔗
2.4.1. Figure (relativement) simple🔗
1#!/usr/bin/env python3
2# Time-stamp: <2023-08-15 12:35:22 ycopin>
3
4"""
5Exemple un peu plus complexe de figure, incluant 2 axes, légendes, axes, etc.
6"""
7
8import numpy as np
9import matplotlib.pyplot as plt
10
11x = np.linspace(-np.pi, 3*np.pi, 2*360)
12
13# Signal carré
14y = np.sign(np.sin(x)) # = ± 1
15
16# 3 premiers termes de la décomposition en série de Fourier
17y1 = 4/np.pi * np.sin(x) # Fondamentale
18y2 = 4/np.pi * np.sin(3*x) / 3 # 1re harmonique
19y3 = 4/np.pi * np.sin(5*x) / 5 # 2de harmonique
20# Somme des 3 premières composantes
21ytot = y1 + y2 + y3
22
23# Figure
24fig = plt.figure() # Création de la Figure
25
26# 1er axe: composantes
27ax1 = fig.add_subplot(2, 1, 1, # 1er axe d'une série de 2 × 1
28 ylabel="Composantes",
29 title="Décomposition en série de Fourier")
30ax1.plot(x, y1, label="Fondamental")
31ax1.plot(x, y2, label="1re harmonique")
32ax1.plot(x, y3, label="2de harmonique")
33ax1.legend(loc="upper left", fontsize="x-small")
34
35# 2nd axe: décomposition
36ax2 = fig.add_subplot(2, 1, 2, # 2d axe d'une série de 2 × 1
37 ylabel="Décomposition",
38 xlabel="x [rad]")
39ax2.plot(x, y, lw=2, color='k', label="Signal carré")
40ax2.plot(x, ytot, lw=2, ls=':', color='k', label="Somme des composantes")
41ax2.legend(loc="upper left", fontsize="x-small")
42
43# Sauvegarde de la figure (pas d'affichage intéractif)
44fig.savefig("figure.png")
Source: figure.py
2.4.2. Filtres du 2nd ordre🔗
Fig. 2.2 Filtre passe-haut du 2nd ordre.🔗
1#!/usr/bin/env python3
2
3import numpy as np
4import matplotlib.pyplot as plt
5
6
7def passeBas(x, Q=1):
8 """
9 Filtre passe-bas en pulsation réduite *x* = omega/omega0, facteur de
10 qualité *Q*.
11 """
12
13 return 1 / (1 - x ** 2 + x / Q * 1j)
14
15
16def passeHaut(x, Q=1):
17
18 return -x ** 2 / (1 - x ** 2 + x / Q * 1j)
19
20
21def passeBande(x, Q=1):
22
23 return 1 / (1 + Q * (x - 1 / x) * 1j)
24
25
26def coupeBande(x, Q=1):
27
28 return (1 - x ** 2) / (1 - x ** 2 + x / Q * 1j)
29
30
31def gainNphase(f, dB=True):
32 """
33 Retourne le gain (éventuellement en dB) et la phase [rad] d'un
34 filtre de fonction de transfert complexe *f*.
35 """
36
37 g = np.abs(f) # Gain
38 if dB: # [dB]
39 g = 20 * np.log10(g)
40 p = np.angle(f) # [rad]
41
42 return g, p
43
44
45def asympGain(x, pentes=(0, -40)):
46
47 lx = np.log10(x)
48 return np.where(lx < 0, pentes[0] * lx, pentes[1] * lx)
49
50
51def asympPhase(x, phases=(0, -np.pi)):
52
53 return np.where(x < 1, phases[0], phases[1])
54
55
56def diagBode(x, filtres, labels,
57 title='', plim=None, gAsymp=None, pAsymp=None):
58 """
59 Trace le diagramme de Bode -- gain [dB] et phase [rad] -- des filtres
60 de fonction de transfert complexe *filtres* en fonction de la pulsation
61 réduite *x*.
62 """
63
64 fig, (axg, axp) = plt.subplots(2, 1, sharex=True) # 2×1 axes
65 axg.set(xscale='log', # Axe des gains
66 ylabel='Gain [dB]')
67 axp.set(xlabel=r'x = $\omega$/$\omega_0$', xscale='log', # Axe des phases
68 ylabel='Phase [rad]')
69
70 lstyles = ['--', '-', '-.', ':']
71 for f, label, ls in zip(filtres, labels, lstyles): # Tracé des courbes
72 g, p = gainNphase(f, dB=True) # Calcul du gain et de la phase
73 axg.plot(x, g, lw=2, ls=ls, label=f"Q={label}") # Gain
74 axp.plot(x, p, lw=2, ls=ls) # Phase
75
76 # Asymptotes
77 if gAsymp is not None: # Gain
78 axg.plot(x, asympGain(x, gAsymp), 'k:', lw=2)
79 if pAsymp is not None: # Phase
80 # axp.plot(x, asympPhase(x,pAsymp), 'k:')
81 pass
82
83 axg.legend(fontsize='small')
84
85 # Labels des phases
86 axp.set_yticks(np.arange(-2, 2.1) * np.pi / 2)
87 axp.set_yticks(np.arange(-4, 4.1) * np.pi / 4, minor=True)
88 axp.set_yticklabels([r'$-\pi$', r'$-\pi/2$', r'$0$', r'$\pi/2$', r'$\pi$'])
89 # Domaine des phases
90 if plim is not None:
91 axp.set_ylim(plim)
92
93 # Ajouter les grilles
94 for ax in (axg, axp):
95 ax.grid() # x et y, majors
96 ax.grid(which='minor') # x et y, minors
97
98 # Ajustements fins
99 gmin, gmax = axg.get_ylim()
100 axg.set_ylim(gmin, max(gmax, 3))
101
102 fig.subplots_adjust(hspace=0.1)
103 axg.xaxis.set_major_formatter(plt.matplotlib.ticker.ScalarFormatter())
104 plt.setp(axg.get_xticklabels(), visible=False)
105
106 if title:
107 axg.set_title(title)
108
109 return fig
110
111
112if __name__ == '__main__':
113
114 x = np.logspace(-1, 1, 1000) # de 0.1 à 10 en 1000 pas
115
116 # Facteurs de qualité
117 qs = [0.25, 1 / np.sqrt(2), 5] # Valeurs numériques
118 labels = [0.25, r'$1/\sqrt{2}$', 5] # Labels
119
120 # Calcul des fonctions de transfert complexes
121 pbs = [ passeBas(x, Q=q) for q in qs ]
122 phs = [ passeHaut(x, Q=q) for q in qs ]
123 pcs = [ passeBande(x, Q=q) for q in qs ]
124 cbs = [ coupeBande(x, Q=q) for q in qs ]
125
126 # Création des 4 diagrammes de Bode
127 figPB = diagBode(x, pbs, labels, title='Filtre passe-bas',
128 plim=(-np.pi, 0),
129 gAsymp=(0, -40), pAsymp=(0, -np.pi))
130 figPH = diagBode(x, phs, labels, title='Filtre passe-haut',
131 plim=(0, np.pi),
132 gAsymp=(40, 0), pAsymp=(np.pi, 0))
133 figPC = diagBode(x, pcs, labels, title='Filtre passe-bande',
134 plim=(-np.pi / 2, np.pi / 2),
135 gAsymp=(20, -20), pAsymp=(np.pi / 2, -np.pi / 2))
136 figCB = diagBode(x, cbs, labels, title='Filtre coupe-bande',
137 plim=(-np.pi / 2, np.pi / 2),
138 gAsymp=(0, 0), pAsymp=(0, 0))
139
140 plt.show()
Source: filtres2ndOrdre.py