Agrégation informatique externe 2023, épreuve 2, option ingénierie informatiqueSujet et rapport du jury

AGREGATION
CONCOURS EXTERNE

MODÉLISATION D'UN SYSTÈME, D'UN PROCÉDÉ OU D'UNE ORGANISATION

INFORMATION AUX CANDIDATS

Ces codes doivent être reportés sur chacune des copies que vous remettrez.

Pilotage d'un exosquelette par une interface Cerveau-Machine

Introduction

Contexte scientifique et technique

Figure 1 - Démonstration de mouvements autonomes contrôlés par la pensée du paraplégique : (a) patient (b) exosquelette (c) installation du patient (d) schéma et radiographie de présentation de l'implant (e) test de marche

CLINATEC

Le projet BCI (Brain Computer Interface)

Système étudié

• -la numérisation des signaux ECoG est assurée par les implants WIMAGINE. Lors d'une opération, deux implants sont posés à demeure dans la boîte crânienne du patient, à la place de l'os crânien et en dessous du cuir chevelu. Ils sont en contact avec le cerveau via des électrodes, pour enregistrer les hémisphères droit et gauche. Ils émettent, via des antennes, un signal reçu par un casque amovible posé sur la tête et connecté à une station de base installée dans le dos de l'exosquelette. Ce signal est synchronisé à la réception par la station de base. Les données sont découpées en lots pour le traitement de l'information en temps réel.

• -Dans un seconde temps, les signaux ECoG sont traités et soumis à des modèles adaptatifs en charge d'estimer l'intention de mouvement du patient. Ce traitement est réalisé sur la cible matérielle BCI PC Board du diagramme de déploiement de la figure 2(b).
• -Les ordres moteurs décodés sont ensuite transmis aux systèmes EMM(EMY Motion Manager) calculant le positionnement recherché pour l'exosquelette et EMC (EMY Motion Controller) gérant les commandes motrices adaptées au pilotage des actionneurs de l'exosquelette.

Figure 2 - (a) Présentation schématisée du projet BCI de Clinatec. Le patient souffrant d'un handicap moteur est placé dans un exosquelette et imagine le mouvement qu'il souhaite exécuter. Un implant cérébral permet de capter l'activité électrique résultante. Puis un système informatique décode ce signal et pilote l'exosquelette. (b) Diagramme de déploiement UML du système BCI.

Figure 3 - Extrapolation du cahier des charges du système BCI développé par CLINATEC

Partie 1. Compression de l'information ECoG pour la transmission implant-station de base

Figure 4 - Système implant et station de base : (a) Chaque implant est constitué d'une matrice de 64 électrodes au-dessus desquelles est positionnée l'électronique embarquée (b) Schéma synthétique de l'implant, contenant notamment un micro-contrôleur MSP420F2618 et deux ASIC CINESIC32 permettant l'enregistrement de 64 canaux en simultané (c) Schéma synthétique des ASIC CINESIC32, on peut noter la présence d'un convertisseur Analogique Numérique (SAR ADC 12 bits) (d) Schéma synthétique de la station de base : le lien UHF (400 MHz) sert à la transmission des données, le lien HF permet la transmission de puissance sans fil vers l'implant.

1.1 Contraintes à la transmission du signal ECoG

• (a)Donner la taille d'un échantillon de tension ECoG pour un seul canal?
• (b)Calculer la taille de l'ensemble d'une mesure des potentiels d'ECoG?

• (c)Calculer le nombre de Data Blocks définis dans la documentation constructeur (document technique DT1) nécessaire pour envoyer l'ensemble de la mesure ECoG.

• (a)Déterminer le débit binaire théorique pour la transmission de l'ensemble des signaux ECoG.
• (b)A partir des données constructeur en document technique DT1, quel est le débit maximum effectif que la liaison UHF peut supporter ? Conclure sur la possibilité de transmettre à priori l'ensemble des données ECoG.
• (c)A f_e = 1kHz, déterminer le nombre de canaux maximum que la liaison sans fil permettra d'envoyer.
• (d)L'ensemble des canaux sont échantillonnés à f_e = 1kHz et transmis. Calculer le nombre de bits maximum sur lesquels doit être codé un échantillon?

• (a)Montrer que 0 ≤ H(s_n^∗).
• (b)A quelle condition sur la séquence s_n^∗ l'entropie H(s_n^∗) est-elle nulle?
• (c)Dans l'hypothèse où chaque caractère de A est équiprobable dans s_n^∗, que vaut l'entropie H(s_n^∗) ?
• (d)Montrer que H(s_n^∗) ≤ N_1.
• (e)Conclure sur le sens physique à donner à la quantité d'entropie H(s_n^∗) de la séquence s_n^∗.

• (a)A partir du prototype donné ci-dessous, écrire le code de la fonction compute_probabilities qui renvoie une variable itérable contenant les probabilités d'apparition des caractères de A dans la séquence s_n^∗. Il est possible d'utiliser la fonction numpy.bincount dont l'aide est donnée en document technique DT3.

def compute_probabilities(sequence, N_bits=12):
    """calcule la probabilite d'apparition de chaque
    mot binaire de N_bits
    Arguments
        sequence : numpy array
            signal ECoG brut sur une electrode en sortie
            du convertisseur Analogique-Numerique
        N_bits : int
            nombre de bits du convertisseur Analogique
    " " "

def compute_entropy(signal, N_bits):
    """calcule l'entropie d'un signal,
    numerise sur un nombre de bits donne
    Arguments
        signal : numpy array
            signal ECoG brut sur une electrode en sortie
            du convertisseurAnalogique-Numerique
        N_bits : int
            nombre de bits du convertisseur Analogique
    " " "

Figure 5 - Histogramme des entropies calculées pour une base de données de signaux ECoG

1.2 Compression spatiale par utilisation de la Transformée en Cosinus Discrète

• (a)Compléter le document réponse DR-1 avec les valeurs numériques des coefficients de C pour N = 8.
• (b)Calculer S(0) en fonction de < s > _(0, N − 1), la moyenne de la séquence s sur N éléments.
• (c)Montrer que
  CC^T = I_N
  où I_N est la matrice identité de taille N. On pourra admettre que :
  • - N est pair,
  • - ∀k ∈ [ [1, N − 1] ], ∑_(i = 0)^(N − 1)cos(((2i + 1)kπ)/(2N)) = 0
• (d)En déduire que C est orthogonale, c'est à dire que C^(− 1) = C^T.
• (e)Retrouver la DCT-1D inverse, c'est-à-dire l'expression de s(k) en fonction de S.
• (f)L'énergie d'une séquence est définie par
  E_S = ‖S‖^2≜∑_(k = 0)^(N − 1)S(k)^2 = S_–^t S_–
  montrer que la DCT-1D conserve l'énergie, c'est-à-dire que E_S = E_s.

Figure 6 - Représentation sous forme de tenseur de dimension 3 des données ECoG, la matrice d'électrodes d'enregistrement est représentée sous forme d'une matrice carrée dont les deux premières dimensions sont spatiales, la troisième dimension représente l'évolution dans le temps

• (a)Implémenter la fonction DCT2D _ matrix prenant en argument les indices i et j et renvoyant la matrice B_(i, j) dont l'entête est donné ci-dessous :

import numpy as np
def DCT2D_matrix(i, j, N=8):
    ,,,
    Construit la matrice B i,j
    Arguments:
    i, j : int
        indices du signal temporel
    Returns:
    B : numpy.array
    ,,,

import numpy as np
def DCT2D_all_matrix(N):
    ,,,
    Construit la liste 2D des matrices B i,j pour une profondeur
    de DCT de N
    Returns:
    mat_listoflist : list of lists
        liste (lignes) de listes (colonnes) des matrices B
        pour une profondeur de DCT de N
    ,,,
    mat_listoflist = []

import numpy as np
def DCT2D(s , k ,N=8):
    ,,,
    Calcule la DCT2D sur l'echantillon temporel
    k du tenseur s
    Arguments :
    s : numpy array
        tenseur ECoG en domaine temporel
    k : int
        indice en temps de l echantillon

N : int
    profondeur de DCT, par defaut a &
Returns:
S : numpy array
    image de s[:,:,k] par la DCT2D
,,,

• (a)En considérant qu'une multiplication coûte 3 périodes d'horloge, calculer le nombre de périodes d'horloge nécessaires pour le calcul de la transformée d'une image spaciale de S.
• (b)Le signal ECoG étant toujours échantillonné à f_e = 1kHz, calculer la fréquence d'horloge minimale pour le micro-contrôleur en ne considérant que le calcul de la DCT.

Figure 7 - Exemple de tenseur du signal ECoG brut et après application de la DCT-2D. (a) tracé de 10 secondes des canaux ECoG (un canal par trace) brut en aval de la conversion analogique numérique (b) même signal, après application de la DCT. Le canal se démarquant correspond au tracé à travers le temps pour les indices (i, j) = (0, 0), l'ensemble des autres indices (i, j) ≠ (0, 0) apparaissent groupés et de plus faibles amplitudes. Pour ce signal et sur la base de données considérée, les valeurs de DCT pour les indices (i, j) = (0, 0) sont bornées entre -4095 et 4096. Les valeurs de DCT pour les indices (i, j) ≠ (0, 0) sont bornées entre -127 et 128.

• (a)Sur combien de bits devrait être codé le canal (i, j) = (0, 0), après l'application de la DCT-2D?
• (b)Sur combien de bits devraient être codés les canaux (i, j) ≠ (0, 0) après l'application de la DCT-2D?
• (c)Calculer le nombre de bits moyens par échantillon.
• (d)Conclure sur la possibilité du système de transmettre l'ensemble des informations ECoG enregistrables dans l'implant (exigence EF4 du diagramme d'exigence), et proposer éventuellement des solutions permettant d'améliorer l'embarquabilité du calcul de la DCT-2D dans l'implant.

Partie 2. Décodage des trajectoires à partir des enregistrements ECoG

Figure 8 - (a) Principe de l'identification d'un décodeur à partir d'enregistrements de signaux ECoG et de mouvements mesurés (b) application du décodeur pour la prédiction d'intention de mouvement.

2.1 Identification des features

• (a)Calculer Ψ(f) la transformée de Fourier de ψ(t).
• (b)Identifier la fréquence centrale de l'ondelette mère.
• (c)Pour un facteur d'échelle s, quelle sera la fréquence centrale de l'ondelette fille centrée en 0(τ = 0) : ψ~((t − τ)/s) = 1/sψ(t/s) ?

• -cette fonction prend en entrée un signal échantillonné sig, la période d'échantillonnage T_e et un vecteur (array_like) contenant les fréquences où l'on cherche à identifier les features,
• -le nombre de valeurs et les valeurs de translations τ sont identiques au nombre d'échantillons et valeurs de temps des échantillons du signal d'entrée,
• -il est possible d'utiliser la fonction morlet, définie précédemment.

import WaveletTransform as wt
import ecog
import numpy as np
# ouverture d un enregistrement ECoG de la base de donnees
# d identification du decodeur (Patient numero 1)
enregistrement = ecog.raw_ECoG('./DB/Patient_1')
T_sample = 1e-3
print(enregistrement.ecog_tensor.shape)
print(type(enregistrement.ecog_tensor[0,0,0]))
# parametres
nsec = 51 # temps de depart
tensor_width = enregistrement.ecog_tensor.shape[0]
start = int ((nsec) / T_sample)
stop = int ((nsec + 1) / T_sample)
freqs = np.linspace(10,150,15) # frequences pour la CWT
# extraction des features entre les secondes 51 et 52
feature_vector = np.zeros((tensor_width**2, stop-start, len(freqs)))
for i in range(enregistrement.tensor_width**2):
    sig = enregistrement.get_electrode_signal(i)[start:stop]
    cwt = np.log10(np.abs(wt.cwt_convolve(sig, freqs, T_sample)))
    feature_vector[i,:,:] = np.transpose(cwt)

(8, 8, 1045828)
<class 'numpy.float64'>

• (a)de la variable sig,
• (b)de la variable cwt,

• (c)de la variable feature_vector

import numpy as np
def cyclic_convolution(f,g):
    N = len ( f)
    conv = np.zeros(N)
    for n in range(N):
        for m in range(N):
            if n-m >= 0:
                conv[n] += f[m]*g[n-m]
            else:
                conv [n] += f [m]*g [N+(n-m)]
    return conv

• (a)
  Calculer la sortie du code suivant :

      import numpy as np
      f = np.array([1,3 ,2])
      g = np.array([3,6,4])
          print(cyclic_convolution(f,g))
  

• (b)Quelle est la complexité algorithmique en temps du calcul de convolution cyclique?
  La Transformée de Fourier Discrète (abrégée DFT) est définie pour un signal échantillonné f par :
  DFT(f)[k] = ∑_(n = 0)^(N − 1)f[n]e^((− 2ȷπ)/Nkn)
  Cette opération est en particulier réalisée par la FFT (Fast Fourier Transform) avec une complexité algorithmique en temps de O(NlogN). La transformée inverse, ou IDFT, est réalisée par la IFFT (Inverse FFT) avec la même complexité.

• (a)Montrer que ∀m ∈ [ [0, N − 1] ] :
  DFT(f_N)[k] = ∑_(n = 0)^(N − 1)f_N[n − m]e^((− 2ȷπ)/Nk(n − m))
• (b)Montrer que
  DFT(f∗_N g)[k] = DFT(f_N)[k] ⋅ DFT(g_N)[k]
• (c)En déduire une méthode de calcul de (f∗_N g).
• (d)Quelle sera la complexité de cette méthode de calcul?

• (a)Proposer une implémentation du calcul de convolution cyclique à partir de l'entête de code suivante :

import numpy.fft.fft as fft
import numpy.fft.ifft as ifft
def FFT_cyclic_convolution(f,g):

import numpy as np
f = np.array([1,3 ,2])
g = np.array([3,6,4])
print(FFT_cyclic_convolution(f,g))

• (a)Écrire le code de la fonction fast_convolution dont l'entête est donné ci-dessous, en réutilisant la fonction de convolution cyclique avec FFT :

import numpy as np
import numpy.fft.fft as fft
import numpy.fft.ifft as ifft
def fast_convolution(f,g):

• (b)Calculer la complexité algorithmique en temps de cette fonction.

import numpy as np
f = np.array([1,3 ,2])
g = np.array([3,6,4])
print(np.convolve(f,g,'same'))
print(np.convolve(f,g,'full'))
print(np.fast_convolution(f,g))

[3 15 28 24 8]
[15 28 24]

• (d)En fonction de N la taille d'origine de f et g comment indexer le résultat de fast_convolution pour obtenir un vecteur de taille N ?

2.2 Décodage des intentions de mouvement

Figure 9 - Mapping spatial, fréquentiel et temporel des données d'enregistrement ECoG pour le décodage des intentions de mouvement [7].

• (a)Calculer la taille d'un vecteur de features x.
• (b)A partir de la taille d'une epoch et le pas δt, calculer le nombre n de vecteurs de features.

• (c)En déduire la taille puis le nombre d'éléments de la matrice X pour la phase de training.
• (d)Les données sont constituées de float64. Calculer le poids mémoire de la matrice X lors de la phase de training.

• (a)Écrire une méthode get_ecog_epoch qui renvoie un tenseur (tableau à 3 dimensions) pour l'epoch commençant au temps ' t '.
• (b)Écrire une méthode compute_features qui pour un temps ' t ' calcule et renvoie le vecteur de feature x. Il est possible d'utiliser la fonction Numpy.ravel comme rappelé dans le document technique DT2 pour transformer le tenseur en vecteur.

• (a)Montrer que
  (∂J)/(∂B) = 2(Y − XB)^T(− X)
• (b)En déduire une solution analytique B^ permettant d'évaluer le modèle au sens des moindres carrés.

Figure 10 - Représentation schématisée des matrices X et Y comme étant la somme de scores (t et u) sur des vecteurs de composantes principales (p et q), ici au nombre de F et de matrices de résidus E et F. La décomposition pour les PLS, contrairement à une décomposition par composantes principales simple a pour but de maximiser la corrélation entre les scores t et u pour dégager des composantes explicatives du modèle φ. Cette figure et les notations sont tirées de [7], attention la notation F désigne un entier, et les vecteurs q sont indicés de 1 à F, F désigne une matrice.

Figure 11 - Pseudo-code et éléments d'interprétation géométrique de l'algorithme NIPALS

• (a)Compléter les calculs des variables dans la boucle d'itération à partir de l'algorithme et du diagramme vu en question précédente. (lignes 15, 17, 21, 24, 27 du DR5)
• (b)Compléter le calcul de la matrice B à la ligne 48 du DR5.

• -un filtre de Kalman,
• -PLS,
• -Multiway (N-way) PLS : NPLS,
• -Penalized regression : SNPLS,
• -Polynomial Penalized PLS : PNPLS.

Partie 3. Commande de l'exosquelette

3.1 Motion manager

3.1.1 Angles d'Euler et problème de singularité.

3.1.2 Quaternions