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, p.ex. `print(a)`
 35        sera résolu comme `a.__str__()`
 36
 37        :return: une chaîne de caractères
 38        """
 39
 40        return f"Animal {'vivant' if self.estVivant else 'mort'}, " \
 41            f"{self.masse:.0f} kg"
 42
 43
 44    def meurt(self):
 45        """
 46        L'animal meurt.
 47        """
 48
 49        self.estVivant = False
 50
 51
 52    def grossit(self, masse):
 53        """
 54        L'animal grossit (ou maigrit) d'une certaine masse (valeur algébrique).
 55
 56        :param float masse: prise (>0) ou perte (<0) de masse.
 57        :raise ValueError: masse non réelle.
 58        """
 59
 60        self.masse += float(masse)
 61
 62
 63# Définition d'une classe héritée ==============================
 64
 65class AnimalFeroce(Animal):
 66    """
 67    Un animal féroce est un animal qui peut dévorer d'autres animaux.
 68
 69    La classe-fille hérite des attributs et méthodes de la
 70    classe-mère, mais peut les surcharger (i.e. en changer la
 71    définition), ou en ajouter de nouveaux:
 72
 73    - la méthode `AnimalFeroce.__init__()` dérive directement de
 74      `Animal.__init__()` (même méthode d'initialisation);
 75    - `AnimalFeroce.__str__()` surcharge `Animal.__str__()`;
 76    - `AnimalFeroce.devorer()` est une nouvelle méthode propre à
 77      `AnimalFeroce`.
 78    """
 79
 80    def __str__(self):
 81        """
 82        Surcharge de la fonction `str()`.
 83        """
 84
 85        return "Animal féroce " \
 86            f"{'bien vivant' if self.estVivant else 'mais mort'}, " \
 87            f"{self.masse:.0f} kg"
 88
 89    def devore(self, other):
 90        """
 91        L'animal (self) devore un autre animal (other).
 92
 93        * Si other est également un animal féroce, il faut que self soit plus
 94          gros que other pour le dévorer. Sinon, other se défend et self meurt.
 95        * Si self dévore other, other meurt, self grossit de la masse de other
 96          (jusqu'à 10% de sa propre masse) et other maigrit d'autant.
 97
 98        :param Animal other: animal à dévorer
 99        :return: prise de masse (0 si self meurt)
100        """
101
102        if isinstance(other, AnimalFeroce) and (other.masse > self.masse):
103            # Pas de chance...
104            self.meurt()
105            prise = 0.
106        else:
107            other.meurt()             # Other meurt
108            prise = min(other.masse, self.masse * 0.1)
109            self.grossit(prise)       # Self grossit
110            other.grossit(-prise)     # Other maigrit
111
112        return prise
113
114
115# Définition d'une autre classe héritée ==============================
116
117class AnimalGentil(Animal):
118    """
119    Un animal gentil est un animal avec un `nom`.
120
121    La classe-fille hérite des attributs et méthodes de la
122    classe-mère, mais peut les surcharger (i.e. en changer la
123    définition), ou en ajouter de nouveaux:
124
125    - la méthode `AnimalGentil.__init__()` surcharge l'initialisation originale
126      `Animal.__init__()`;
127    - `AnimalGentil.__str__()` surcharge `Animal.__str__()`;
128    """
129
130    def __init__(self, masse, nom='Youki'):
131        """
132        Initialisation d'un animal gentil, avec son masse et son nom.
133        """
134
135        # Initialisation de la classe parente (nécessaire pour assurer
136        # l'héritage)
137        Animal.__init__(self, masse)
138
139        # Attributs propres à la classe AnimalGentil
140        self.nom = nom
141
142    def __str__(self):
143        """
144        Surcharge de la fonction `str()`.
145        """
146
147        return f"{self.nom}, un animal gentil " \
148            f"{'bien vivant' if self.estVivant else 'mais mort'}, " \
149            f"{self.masse:.0f} kg"
150
151    def meurt(self):
152        """
153        L'animal gentil meurt, avec un éloge funéraire.
154        """
155
156        Animal.meurt(self)
157        print(f"Pauvre {self.nom} meurt, paix à son âme...")
158
159
160if __name__ == '__main__':
161
162    # Exemple d'utilisation des classes définies ci-dessus
163
164    print("Une tragédie en trois actes".center(70, '='))
165
166    print("Acte I: la vache prend 10 kg.".center(70, '-'))
167    vache = Animal(500.)        # Instantiation d'un animal de 500 kg
168    vache.grossit(10)           # La vache grossit de 10 kg
169    print(vache)
170
171    print("Acte II: Dumbo l'éléphant".center(70, '-'))
172    elephant = AnimalGentil(1000., "Dumbo")  # Instantiation d'un animal gentil
173    print(elephant)
174
175    print("Acte III: le féroce lion".center(70, '-'))
176    lion = AnimalFeroce(200)    # Instantiation d'un animal féroce
177    print(lion)
178
179    print("Scène tragique: le lion dévore l'éléphant...".center(70, '-'))
180    lion.devore(elephant)       # Le lion dévore l'éléphant
181
182    print(elephant)
183    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        import sys
 85
 86        eps = sys.float_info.epsilon  # Le plus petit float non nul
 87
 88        return (abs(self.x) <= eps) and (abs(self.y) <= eps)
 89
 90    def distance(self, other):
 91        """
 92        Méthode de calcul de la distance du point (`self`) à un autre point
 93        (`other`).
 94
 95        >>> A = Point(1,0); B = Point(1,1); A.distance(B)
 96        1.0
 97        """
 98
 99        # hypot(dx, dy) = sqrt(dx**2 + dy**2)
100        return ((self.x - other.x) ** 2 + (self.y - other.y) ** 2) ** 0.5
101
102
103# Définition du point origine O
104O = Point(0, 0)
105
106
107# Héritage de classe ==============================
108
109
110class Vector(Point):
111
112    """
113    Un `Vector` hérite de `Point` avec des méthodes additionnelles
114    (p.ex. la négation d'un vecteur, l'addition de deux vecteurs, ou
115    la rotation d'un vecteur).
116    """
117
118    def __init__(self, A, B):
119        """
120        Définit le vecteur `AB` à partir des deux points `A` et `B`.
121
122        >>> Vector(Point(1, 0), Point(1, 1))  # doctest: +ELLIPSIS
123        <....Vector object at 0x...>
124        >>> Vector(0, 1)
125        Traceback (most recent call last):
126        ...
127        AttributeError: 'int' object has no attribute 'x'
128        """
129
130        # Initialisation de la classe parente
131        Point.__init__(self, B.x - A.x, B.y - A.y)
132
133        # Attribut propre à la classe dérivée
134        self.sqnorm = self.x**2 + self.y**2  # Norme du vecteur au carré
135
136    def __str__(self):
137        """
138        Surcharge de la fonction `str()`: `str(self)` sera résolu comme
139        `Vector.__str__(self)` (et non pas comme `Point.__str__(self)`)
140
141        >>> str(Vector(O, Point(0, 1)))
142        'Vector (x=0.0, y=1.0)'
143        """
144
145        return f"Vector (x={self.x}, y={self.y})"
146
147    def __eq__(self, other):
148        """
149        Surcharge du test d'égalité `{self}=={other}`: l'instruction sera
150        résolue comme `self.__eq__(other)`.
151
152        >>> Vector(O, Point(0, 1)) == Vector(Point(1, 0), Point(1, 1))
153        True
154        """
155
156        # On teste ici la nullité de la différence des 2 vecteurs.
157        # D'autres tests auraient été possibles -- égalité des
158        # coordonnées, nullité de la norme de la différence, etc. --
159        # mais on tire profit de la méthode héritée `Point.is_origin()`
160        # testant la nullité des coordonnées (à la précision numérique
161        # près).
162        return (self - other).is_origin()
163
164    def __add__(self, other):
165        """
166        Surcharge de l'opérateur binaire `{self} + {other}`: l'instruction
167        sera résolue comme `self.__add__(other)`.
168
169        On construit une nouvelle instance de `Vector` à partir des
170        coordonnées propres à l'objet `self` et à l'autre opérande
171        `other`.
172
173        >>> A = Point(1, 0); B = Point(1, 1)
174        >>> Vector(A, B) + Vector(B, O) == Vector(A, O)
175        True
176        """
177
178        return Vector(O, Point(self.x + other.x, self.y + other.y))
179
180    def __sub__(self, other):
181        """
182        Surcharge de l'opérateur binaire `{self} - {other}`: l'instruction
183        sera résolue comme `self.__sub__(other)`.
184
185        Attention: ne surcharge pas l'opérateur unaire `-{self}`, géré
186        par `__neg__`.
187
188        >>> A = Point(1, 0); B = Point(1, 1)
189        >>> str(Vector(A, B) - Vector(A, B))  # Différence
190        'Vector (x=0.0, y=0.0)'
191        >>> -Vector(A, B)                     # Négation
192        Traceback (most recent call last):
193        ...
194        TypeError: bad operand type for unary -: 'Vector'
195        """
196
197        return Vector(O, Point(self.x - other.x, self.y - other.y))
198
199    def __abs__(self):
200        """
201        Surcharge la fonction `abs()` pour retourner la norme du vecteur.
202
203        >>> abs(Vector(O, Point(1, 0)))
204        1.0
205        """
206
207        # On pourrait utiliser sqrt(self.sqnorm), mais c'est pour
208        # illustrer l'utilisation de la méthode héritée
209        # `Point.distance`...
210        return Point.distance(self, O)
211
212    def rotate(self, angle, deg=False):
213        """
214        Rotation (dans le sens trigonométrique) du vecteur par un `angle`,
215        exprimé en radians ou en degrés.
216
217        >>> Vector(O, Point(0, 1)).rotate(90, deg=True) == Vector(O, Point(-1, 0))
218        True
219        """
220
221        from cmath import rect  # Bibliothèque de fonctions complexes
222
223        # On calcule la rotation en passant dans le plan complexe
224        z = complex(self.x, self.y)
225        phase = angle if not deg else angle / 57.29577951308232  # [rad]
226        zu = z * rect(1.0, phase)  # Rotation complexe (* exp(i*phase))
227
228        return Vector(O, Point(zu.real, zu.imag))
229
230
231def circumscribedCircle(M, N, P):
232    """
233    Calcule le centre et le rayon du cercle circonscrit aux points
234    M, N, P.
235
236    Retourne: (centre [Point], rayon [float])
237
238    Lève une exception `ValueError` si le rayon ou le centre du cercle
239    circonscrit n'est pas défini.
240
241    >>> M = Point(-1, 0); N = Point(1, 0); P = Point(0, 1)
242    >>> C, r = circumscribedCircle(M, N, P)  # Centre O, rayon 1
243    >>> C.distance(O), round(r, 6)
244    (0.0, 1.0)
245    >>> circumscribedCircle(M, O, N)         # Indéfini
246    Traceback (most recent call last):
247    ...
248    ValueError: Undefined circumscribed circle radius.
249    """
250
251    MN = Vector(M, N)
252    NP = Vector(N, P)
253    PM = Vector(P, M)
254
255    # Rayon du cercle circonscrit
256    m = abs(NP)  # |NP|
257    n = abs(PM)  # |PM|
258    p = abs(MN)  # |MN|
259
260    d = (m + n + p) * (-m + n + p) * (m - n + p) * (m + n - p)
261    if d > 0:
262        rad = m * n * p / d**0.5
263    else:
264        raise ValueError("Undefined circumscribed circle radius.")
265
266    # Centre du cercle circonscrit
267    d = -2 * (M.x * NP.y + N.x * PM.y + P.x * MN.y)
268    if d == 0:
269        raise ValueError("Undefined circumscribed circle center.")
270
271    om2 = Vector(O, M).sqnorm  # |OM|**2
272    on2 = Vector(O, N).sqnorm  # |ON|**2
273    op2 = Vector(O, P).sqnorm  # |OP|**2
274
275    x0 = -(om2 * NP.y + on2 * PM.y + op2 * MN.y) / d
276    y0 = (om2 * NP.x + on2 * PM.x + op2 * MN.x) / d
277
278    return (Point(x0, y0), rad)  # (centre [Point], R [float])
279
280
281if __name__ == "__main__":
282
283    # start-argparse
284    import argparse
285
286    parser = argparse.ArgumentParser(
287        usage="%(prog)s [-p/--plot] [-i/--input coordfile | x1,y1 x2,y2 x3,y3]",
288        description="Compute the circumscribed circle to 3 points in the plan.",
289    )
290    parser.add_argument(
291        "coords",
292        nargs="*", type=str, metavar="x,y",
293        help="Coordinates of point"
294    )
295    parser.add_argument(
296        "-i", "--input",
297        nargs="?", type=argparse.FileType("r"),
298        help="Coordinate file (one 'x,y' per line)",
299    )
300    parser.add_argument(
301        "-P", "--plot",
302        action="store_true", default=False,
303        help="Draw the circumscribed circle",
304    )
305    parser.add_argument(
306        "-T", "--tests",
307        action="store_true", default=False,
308        help="Run doc tests"
309    )
310    parser.add_argument("--version", action="version", version=__version__)
311
312    args = parser.parse_args()
313    # end-argparse
314
315    if args.tests:  # Auto-test mode
316        import sys
317        import doctest
318
319        fails, tests = doctest.testmod(verbose=True)  # Run doc tests
320        sys.exit(fails > 0)
321
322    if args.input:  # Lecture des coordonnées du fichier d'entrée
323        # Le fichier a déjà été ouvert en lecture par argparse (type=file)
324        args.coords = [coords
325                       for coords in args.input
326                       if not coords.strip().startswith("#")]
327
328    if len(args.coords) != 3:  # Vérifie le nb de points
329        parser.error("Specify 3 points by their coordinates 'x,y' "
330                     f"(got {len(args.coords)})")
331
332    points = []                 # Liste des points
333    for i, arg in enumerate(args.coords, start=1):
334        try:                    # Déchiffrage de l'argument 'x,y'
335            x, y = (float(t) for t in arg.split(","))
336        except ValueError:
337            parser.error(f"Cannot decipher coordinates #{i}: {arg!r}")
338
339        points.append(Point(x, y))      # Création du point et ajout à la liste
340        print(f"#{i:d}: {points[-1]}")  # Affichage du dernier point
341
342    # Calcul du cercle cisconscrit (lève une ValueError en cas de problème)
343    center, radius = circumscribedCircle(*points)  # Délistage
344    print(f"Circumscribed circle: {center}, radius: {radius}")
345
346    if args.plot:  # Figure
347        import matplotlib.pyplot as plt
348
349        fig, ax = plt.subplots()
350        ax.set(aspect="equal")
351        # Points
352        ax.plot([p.x for p in points], [p.y for p in points], "ko")
353        for i, p in enumerate(points, start=1):
354            ax.annotate(f"#{i}", (p.x, p.y),
355                        xytext=(5, 5), textcoords="offset points")
356        # Cercle circonscrit
357        c = plt.matplotlib.patches.Circle((center.x, center.y),
358                                          radius=radius, fc="none", ec="k")
359        ax.add_patch(c)  # Cercle
360        ax.plot(center.x, center.y, "r+")  # Centre
361
362        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
Cercle circonscrit.

Source: circonscrit.py

2.4. Matplotlib🔗

2.4.1. Figure (relativement) simple🔗

Exemple de figure simple.

Figure 2.1 Exemple de figure (charte graphique: seaborn)🔗

 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🔗

Filtre passe-haut du 2nd ordre

Figure 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