[PYTHON] Une introduction à l'analyse vocale pour les applications musicales

Ceci est un résumé du contenu de Traitement du signal audio pour les applications musicales que j'ai pris à Coursera. Je pense que ce sera une passerelle lorsque vous voudrez aller plus loin que la simple lecture d'une piste.

Objectif de l'analyse vocale

En gros, l'analyse vocale consiste à clarifier "de quelle fréquence le son que vous entendez est composé". C'est ce qu'on appelle l'analyse spectrale. Une fois que cela sera révélé, vous pourrez:

Caractériser à partir de la composition du son (identifier l'instrument / le haut-parleur, estimer le genre de musique, etc.)
Exclure les sons autres que la fréquence souhaitée (filtrage)
Composez le son souhaité en synthétisant des sons (synthétiseur, etc.)

Comme mentionné ci-dessus, si la décomposition (analyse) est réussie, il sera possible de créer un son en synthétisant le contraire. En d'autres termes, l'analyse musicale est un processus de ** décomposition-> analyse-> (conversion / filtrage) -> reconstruction ** dans son ensemble (on peut dire qu'il s'agit d'un alchimiste d'acier).

Expression de la voix

La première étape de l'analyse de la parole commence par l'expression de la parole en tant que fonction.

x[n] = Acos(\omega nT + \phi) = Acos(2\pi fnT + \phi)

$ A $: amplitude amplitude-> volume
$ \ omega $: fréquence angulaire dans le rayon Fréquence angulaire-> Pas (fréquence)
f: frequency in Hz (\frac{w}{2\pi})
$ \ phi $: phase initiale dans le rayon Phase initiale
n: time index
$ f_s $: fréquence d'échantillonnage du taux d'échantillonnage
T: sampling period in second (1/f_s, t=nT=n/f_s)

Je pense que c'est dans ma mémoire que le graphique de la fonction triangulaire ressemble à une onde sonore. Ceci est utilisé pour exprimer le son en tant que fonction. La vitesse angulaire de $ \ omega $, comme son nom l'indique, est l'angle auquel elle se déplace par unité de temps (1 seconde). L'unité d'angle est le radian (pour ceux qui disent radian, veuillez vous reporter à ici), dans le cas du radian $ 2 \ pi $ fait un tour, c'est-à-dire un cycle, donc en divisant $ \ omega $, dont l'unité est le radian, par $ 2 \ pi $ donne $ f $ par seconde, soit Hz.

Le point à noter est $ n , qui est l'indice du son final ( x [n] $). Si ce $ n $ avance de 1, il avance de $ T $ secondes. $ T $ est l'inverse de la fréquence d'échantillonnage, qui est l'intervalle d'échantillonnage. La fréquence d'échantillonnage représente la fréquence avec laquelle le son est acquis, et plus elle est élevée, plus le son (son de fréquence plus élevée) peut être traité. La fréquence d'échantillonnage doit avoir une marge d'au moins deux fois la fréquence que vous souhaitez acquérir pour la composition (théorème d'échantillonnage), et même si la plage audible humaine des CD généraux est d'environ 20 000 Hz. 44 100 Hz, soit environ le double, est adopté.

La fréquence d'échantillonnage de la source sonore appelée haute résolution est de 48 000 Hz, ce qui est plus élevé que le taux des CD généraux, raison pour laquelle elle est appelée haute résolution = haute résolution. Si vous ne pouvez pas l'entendre, vous pourriez penser que cela n'a pas de sens, mais il semble que cela inclut la sensation d'air que vous ressentez avec votre corps au lieu de vos oreilles.

L'introduction a été allongée, mais la figure ci-dessous est un graphique de cette fonction avec les paramètres réels. Lorsque vous l'entendez comme un son, cela produit un son semblable à celui d'un test auditif.

CC by MTG(Basicmathematics-Sinewaveplot)

Cependant, il est difficile à gérer tel quel, nous allons donc utiliser la formule d'Euler pour le convertir en fonction exponentielle.

\bar{x}[n] = Ae^{j(\omega nT + \phi)} = Acos(\omega nT + \phi) + Asin(\omega nT + \phi)

La formule d'Euler montre la relation entre l'exposant complexe et la fonction triangulaire $ e ^ {j \ phi} = cos \ phi + jsin \ phi $ (voir la figure ci-dessous), et en appliquant cela, la fonction triangulaire Il est possible de convertir la représentation en une représentation d'un exposant complexe.

CC by MTG(Basicmathematics-Euler'sformula)

L'avantage d'utiliser la représentation exponentielle complexe est que l'amplitude $ A $ peut être facilement exprimée comme la valeur absolue du nombre complexe, et la phase initiale $ \ phi $ peut être facilement exprimée comme l'angle de déviation.

CC by MTG(Basicmathematics-Complexnumbers)

L'histoire jusqu'à présent peut être résumée comme suit.

Supposons que le son est une fonction périodique liée au temps et exprimez-le avec une fonction triangulaire.
Plus l'amplitude $ A $ est grande, plus le son est fort et plus la vitesse angulaire $ \ omega $ est élevée, plus le son est élevé.
Pour faciliter l'utilisation, la représentation de fonction triangulaire est convertie en une représentation d'exposant complexe en utilisant la formule d'Euler.
Avec l'introduction de la représentation exponentielle complexe, l'amplitude peut être représentée par une valeur absolue et la phase initiale peut être représentée par un angle de déviation.
Le son est exprimé sous la forme d'une valeur par unité de temps déterminée par la fréquence d'échantillonnage (expression de temps discret). Plus précisément, l'intervalle de temps entre $ x [n] $ et $ x [n + 1] $ est $ 1 / f_s $ lorsque la fréquence d'échantillonnage est $ f_s $.

Maintenant que le son peut être exprimé, je vais vous expliquer comment décomposer et reconstruire ce son.

Décomposition / reconstruction de la voix

Transformée de Fourier discrète (DFT)

Dès la conclusion, c'est la transformée de Fourier qui décompose la voix et la transformée de Fourier inverse qui la reconstruit. Cependant, comme le temps géré par l'ordinateur est discret comme expliqué ci-dessus, il est appelé transformation de Fourier discrète / transformation inverse.

Tout d'abord, voyons comment la transformation de Fourier discrète (DFT) qui effectue la décomposition est définie.

X[k] = \sum^{N-1}_{n=0} x[n] e^{-j2\pi kn/N}

$ N $: Nombre d'échantillons
n: discrete time index
k: discrete frequency index k=0,...,N-1
$ \ omega_k = 2 \ pik / N $: fréquence angulaire
$ f_k = f_s k / N $: Fréquence (Hz)

L'essence de la transformation de Fourier est de découvrir la répartition $ X [k] $ de la quantité de son que chaque fréquence $ k $ est contenue dans les échantillons $ N $ obtenus $ x [n] $. Cependant, il faut noter ici que $ k $ n'est pas une simple fréquence. La limite supérieure de la fréquence qui peut être examinée est déterminée par la fréquence d'échantillonnage $ f_s $, mais si seulement $ N $ est échantillonné à partir de celle-ci, elle sera examinée tous les $ f_s / N $ (= discret). Ce $ f_s / N $ est l'unité de k, donc la relation entre $ k $ et la fréquence est $ f_k = f_s k / N $ comme décrit ci-dessus. Le fait est que la finesse de la fréquence à examiner dépend du nombre d'échantillons.

$ n $ est une unité de temps discret, $ 1 / f_s $ secondes
$ k $ est une unité de fréquence discrète, unité $ f_s / N $ Hz (cela s'appelle bin)

Si vous ne gardez pas ces deux points à l'esprit, vous serez confus par la conversion d'unité. La partie $ e ^ {-j2 \ pi kn / N} $ dans DFT a pour effet d'annuler tout autre chose que la fréquence souhaitée $ k $ à partir de $ x [n] $ (The Discrete Fourier Transform 2 of 2, [Généralisations d'identité d'Euler](http://en.wikipedia.org/ Voir wiki / Euler's_identity)), vous ne pouvez obtenir que $ k $ en multipliant $ x [n] $ par ceci. En calculant ceci pour $ n $ à chaque fois et en prenant la somme, nous calculons la somme des composants $ k $ contenus dans $ n $ à chaque fois, c'est-à-dire combien de fréquence $ k $ est incluse dans tout le son. Tu peux le faire.

La figure ci-dessous montre ce qui a été réellement calculé et tracé.

CC by MTG(TheDiscreteFourierTransform1of2)

La première étape montre le son réel, les deuxième et troisième étapes montrent le résultat DFT, la deuxième étape est l'amplitude et la troisième étape est le graphique de phase initiale. Dans la deuxième rangée, l'axe vertical est l'intensité de l'amplitude convertie en digibel, et l'axe horizontal est Hz, qui montre la hauteur du son et sa force (c'est ce qu'on appelle le spectre de magnitude). ). Dans la troisième rangée, l'axe vertical est la phase initiale (vitesse angulaire) et l'axe horizontal est Hz, ce qui indique quel son de hauteur commence à sonner à quel moment (c'est ce qu'on appelle le spectre de phase).

Ici, je résumerai quelques questions connexes.

Symétrie En analysant le son réel, la magnitude / phase est symétrique. C'est $ A_ocos (2 \ pi k_0 n / N) = \ frac {A_0} {2} e ^ {j2 \ pi k_0 n / M} + \ frac {A_0} {2} e ^ {-j2 \ pi k_0 Il peut être dérivé en développant la formule de n / M} $, et dans ce cas, la grandeur apparaît sous la forme d'une symétrie de ligne centrée sur les deux points de $ k = k_0 et -k_0 $. (Strictement parlant, dans le cas d'un son composé uniquement de la partie réelle (son ordinaire), la magnitude est un objet pair et la phase est étrangement symétrique, et lorsqu'elle est composée de la partie réelle et de la partie imaginaire, la magnitude est un objet pair et la phase est 0)

CC by MTG(TheDiscreteFourierTransform2of2DFTofrealsinusoid)

Puisque nous n'avons pas besoin du côté négatif, nous ne traitons généralement que 0 ou plusieurs parties (moitié positive).

Energy
La somme des carrés de grandeur est appelée énergie.

\sum_{n=-N/2}^{N/2-1}|x[n]|^2 = \frac{1}{N}\sum_{k=-N/2}^{N/2-1}|X[k]|^2

Décibel (décibel) Une unité qui représente la magnitude et est définie ci-dessous.

db = 20 * log_{10}(abs(X)=magnitude)

Phase unwrapping
La phase est $ 2 \ pi $, c'est-à-dire qu'elle revient à 0 après une rotation, donc le traitement de la phase avec cet ajouté s'appelle le déballage. En effectuant ce processus, il devient plus facile de voir la transition de la phase.

CC by MTG(FourierTransformproperty2of2Phaseunwrapping)

Zero-padding
En ajoutant un élément de voix 0 après la voix à analyser, la courbe d'amplitude obtenue par DFT devient lisse et le pic peut être facilement estimé. Par conséquent, le nombre d'échantillons $ N $, qui est la taille de DFT, est souvent plus grand que la taille réelle du son $ M $, et la différence est souvent remplie de 0.

CC by MTG(FourierTransformproperty2of2Zero-padding)

Effectuer une DFT et clarifier la magnitude et la phase est le processus équivalent à la «décomposition sonore».

Reconstruction (IDFT: Inverse Discrete Fourier Transform)

La transformée de Fourier discrète inverse (IDFT) est utilisée pour reconstruire le son original en utilisant la magnitude et la phase obtenues à la suite de la DFT.

x[n] = \frac{1}{N}\sum^{N-1}_{k=0}X[k]e^{j2\pi kn/N}

Il s'agit de l'opération inverse de DFT et resynthétise la voix.

Accélérant

Lors du calcul, utilisez la transformation de Fourier à grande vitesse / transformation inverse (FFT / IFFT) qui accélère la DFT / IDFT. Pour utiliser FFT, procédez comme suit:

Fenêtre à phase nulle (voir la figure ci-dessous car il est difficile à expliquer avec des mots)
Faites correspondre la taille à la puissance de 2 (128, 256, 512 ...)
Remplissez la mauvaise taille avec 0 (Zero-padding)

FFT / IFFT est implémenté dans scipy et peut être facilement implémenté en l'utilisant.

L'explication ci-dessus est plus rapide à lire dans le code, veuillez donc vous y référer également.

MTG/sms-tools/software/models/dftModel.py

Je pense que vous pouvez maintenant décomposer et reconstruire la voix. Cependant, dans l'état actuel des choses, nous reconstruisons simplement le son original, je voudrais donc aller un peu plus loin sur la façon d'analyser les résultats de la décomposition à partir d'ici, puis enfin regarder les parties liées à la conversion et au filtrage.

Analyse des résultats de la décomposition vocale

Analyse des changements de séries chronologiques (STFT: Short-Time Fourier Transform)

Une chanson typique a une durée de plusieurs minutes, et lorsqu'elle est appliquée à la FFT, les sons qui sonnent à différents moments de la chanson sont mélangés, et l'ampleur et la phase deviennent encombrées et les caractéristiques deviennent difficiles à saisir. Je vais finir. Par conséquent, le morceau est divisé en longueurs fixes et une transformée de Fourier discrète est exécutée pour chaque unité de division. En les organisant, vous pouvez saisir le changement de magnitude / phase dans les séries chronologiques. Cette méthode est appelée transformation de Fourier à court terme (STFT).

X_l[k] = \sum^{N/2-1}_{n=-N/2}w[n]x[n+lH]e^{-j2\pi kn/N}

w: analysis window
l: frame number (0,1,...)
H: hop size

Il s'agit d'une image de décalage de la plage d'analyse de taille N par H (voir la figure ci-dessous).

CC by MTG(TheShort-TimeFourierTransform(1of2))

$ w $ est appelé une fonction de fenêtre. La raison de l'application d'une telle chose est qu'on suppose que la partie découpée de taille N est un "signal périodique". Cependant, en réalité ce n'est pas le cas, donc la première partie est appliquée avec une fonction qui converge vers 0 à la fin pour faire apparaître comme s'il s'agissait d'un signal périodique. La fonction pour cela s'appelle la fonction de fenêtre.

Vous vous demandez peut-être s'il est correct d'appliquer quelque chose comme ça, mais au contraire, si vous n'appliquez pas la fonction de fenêtre, l'hypothèse de "signal périodique" sera rompue, donc lorsque vous reconstruisez la voix, du bruit sera généré. Il semblera qu'il est inclus.

Plus précisément, la fonction de fenêtre est la suivante (ce qui suit est une simple fenêtre rectangulaire).

CC by MTG(TheShort-TimeFourierTransform(1of2)Analysiswindow)

Il prend une valeur de 1 au pic et prend la forme d'un seul pic qui se désintègre autour de lui. En prenant la magnitude, il devient comme indiqué dans la figure ci-dessus. Le pic est appelé lobe principal et le pic à côté du pic est appelé lobe latéral. Puisqu'il est préférable que la fonction de fenêtre passe uniquement la fréquence cible, plus le lobe principal est étroit, plus la résolution de fréquence est élevée et plus le lobe latéral autre que le lobe principal est bas, plus l'effet sur les sons (petits sons) autres que le principal est important. Va être réduit. Il y a un compromis entre l'étroitesse du lobe principal et du lobe latéral bas, et il est nécessaire de les utiliser correctement en fonction de la situation. En général, sélectionnez celui avec un lobe latéral bas dans le cas où des sons forts et faibles sont inclus, et sélectionnez celui avec un lobe principal étroit dans le cas où il tient dans un certain volume, en donnant la priorité à la résolution. (Le premier est souvent blackman / blackman-harris et le second martèle / pend).

Référence
Window function　Windowing Understanding FFT Windows Choosing a Windowing Function

Sinusoidal model

Le modèle sinusoïdal consiste à penser un son complexe comme un ensemble de sons simples (sinusoïde) d'une fréquence spécifique. Dans le modèle ci-dessous, le son au temps n est représenté par R sons de fréquences spécifiques.

y[n] = \sum^R_{r=1} A_r[n]cos(2\pi f_r [n]n)

$ R $: Nombre d'ondes sinusoïdales
$ A_r [n] $: Fluctuation au temps n
$ f_r [n] $: Fréquence au temps n

S'il peut s'appliquer à ce modèle, il semble possible de reconstruire des sons complexes à partir de sons simples. Alors, comment savoir de quelle fréquence est composé un son? L'indice est toujours le spectre.

L'analyse spectrale révèle que les sons de différentes fréquences sont détectés sous forme de pics de magnitude, comme indiqué ci-dessous.

Week 5: Sinusoidal model, Theory lecture 1: Sinusoidal model 1

Par conséquent, il semble que les fréquences qui composent le son peuvent être bien identifiées en suivant les étapes ci-dessous.

Créez un spectre de magnitude afin que les pics puissent être bien détectés
Identifiez l'emplacement du pic

Commencez par aller de "1. Créez un spectre de magnitude pour que les pics puissent être bien détectés". Le point de ceci est la taille de la fenêtre, et il est nécessaire de définir de sorte que plusieurs composants de fréquence tiennent à l'intérieur de la fenêtre.

Week 5: Sinusoidal model, Theory lecture 1: Sinusoidal model 1

Comme indiqué dans ci-dessus, la fréquence détectable est en unités $ f_s / N . Par conséquent, cela devrait être inférieur à la différence entre les deux fréquences que vous essayez de détecter (car sinon vous ne pouvez pas décomposer les deux fréquences). Tout d'abord\frac{f_s}{|f_{k+1} - f_k|}Est dérivé. Et le nombre de bacs qui s'insèrent dans le lobe principal de la fonction de fenêtre est fixe, et c'estB_sSera. Donc,B_s \frac{f_s}{|f_{k+1} - f_k|}$La taille appropriée peut être dérivée avec.

Un exemple d'un mélange des 440/490 ci-dessus est le suivant.

Week 5: Sinusoidal model, Theory lecture 1: Sinusoidal model 1

Eh bien, je ne sais pas vraiment quelle est la fréquence de ce son (440 et 490!), Donc je vais ajouter quelque chose comme 100-2000Hz et choisir un M qui fonctionne bien sur toute la fréquence entre les deux. (Comme ça. L'axe horizontal est la fréquence et l'axe vertical est k (M = k * 100 +) Puisque k = 21 est stable à n'importe quelle fréquence, 2101 est utilisé comme M).

Comment détecter les pics
Décomposition / synthèse par modèle sinusoïde

Harmonic model Sinusoidal plus residual model(Stochastic model)

Conversion / filtrage

filtration

**A3: Fourier Properties Part-4: Suppressing frequency components using Citation du modèle DFT **

Implementation

conversion

Week8 - Sound transformations

Classification des voix

Week9 - Sound and music description

Vers le monde de la recherche en analyse vocale

Week10 - Concluding topics