Acoustic passion

Introduction

J'ai mis au point un programme simple sous Octave permettant de générer un filtre de correction FIR pour effectuer une correction.

J'effectue pour ma part une correction en champ proche (~30-40cm avec un micro EMM6 avec incidence rasante, c'est à dire 30-45° environ par rapport à la source) et une correction supplémentaire au point d'écoute avec un égaliseur IIR pour compenser les effets de la pièce : c'est la méthode que je préconise.

Installation

Tout d'abord il faut installer Octave, voici un lien vers un tutoriel pour l'installation d'Octave sur windows : http://enacit1.epfl.ch/cours_matlab/octave.shtml#windows

Vous pouvez suivre la "Procédure d'installation rapide" pour installer simplement et rapidement Octave (il suffit de dézipper le contenu dans la racine C:/ et déplacer un raccourci sur le bureau... j'espère ne pas vous perdre ici !).

A présent vous pouvez télécharger mon logiciel ici : https://drive.google.com/file/d/0B_qpfrOOECB0TnFDWURYNkpvRG8/edit?usp=sharing

Dézipper le contenu quelque part, lancez Octave puis rendez vous dans la racine de mon logiciel. Pour cela vous devez entrer la commande "cd " puis le chemin vers le répertoire du logiciel. Si vous l'avez dézipper dans C:/ par exemple, vous devez taper "cd C:/Octave-DSP" puis valider par entrée. Vous pouvez copier/coller le chemin du répertoire ou encore glisser/déposer le répertoire dans la fenêtre d'octave pour entrer le chemin d'accès immédiatement.

Vous devriez être dans le répertoire du logiciel, vérifiez en tapant "ls" puis entrée. Si vous êtes dans le répertoire du programme vous devriez voir apparaître la liste des fichiers qui sont contenus dans ce répertoire.

Les fichiers / fonctions :

Le programme est scindé en plusieurs fonctions réparties dans les différents fichiers. Voici la fonctionnalité de chaque fonction :

• brickwall_filter.m => function filter=brickwall_filter(length, highpass, lowpass)

Fonction qui prend une longueur (en samples), une fréquence de coupure passe-haut et passe-bas et génère un filtre brickwall avec fenêtrage blackman.

• dBToLinear.m => function result=dbToLinear(value)

Fonction qui prend une valeur en dB et la converti en linéaire (avec la relation dB=20*log10(linéaire))

• detect_rejectband.m => function [highpass lowpass]=detect_rejectband(fftSignal, threshold)

Fonction qui prend un signal sous forme FFT, et une valeur de dynamique en dB (exemple : -80dB) et qui renvoie les fréquences autour du pic maximum à partir de laquelle l'amplitude descend en dessous de cette dynamique. Par exemple si le signal est un filtre passe-bande de -5dB et qu'on utilise une dynamique de -80dB, le programme partira d'une valeur maximale (donc dans la bande de passage du filtre), et ira cherchez du côté gauche la fréquence de passe-haut à partir de laquelle l'amplitude < -85dB. De même pour le côté droit. La fonction renvoie alors les 2 fréquences passe-haut et passe-bas.

• fftdiv.m => function result=fftdiv(signal1, signal2)

Fonction qui prend 2 signaux (domaine temporel) et qui calcule le signal "result" tel que result * signal2 = signal1. Comme dans la convolution A*B=C il faut que la longueur de C=longueur de A+B-1, la taille du signal "result" obtenu vaut longueur signal1-signal2+1. Il faut donc que signal1 > signal2. Evidemment plus le signal "result" calculé est grand, meilleure est la correction.

• generate_correction.m => function correction=generate_correction(measure, target)

Fonction qui prend un signal de mesure (la réponse d'impulsion du système en fait, dans le domaine temporel donc) et un signal cible (la réponse d'impulsion souhaitée, dans le domaine temporel encore) et qui renvoie la réponse d'impulsion de la correction à convolutionner au système pour atteindre la réponse cible souhaitée.

• max_absolute_peak.m => function index=max_absolute_peak(signal)

Fonction qui prend un signal et qui renvoie l'index du pic maximal en valeur absolue.

• normalize_freq_peak => function result=normalize_freq_peak(signal)

Fonction qui prend un signal (domaine temporel), calcule la transformée de Fourrier et détermine la valeur maximale d'amplitude dans le domaine fréquentiel. Le signal est ensuite amplifié ou atténué de sorte que le pic maximal dans le domaine fréquentiel corresponde à 0dBFS (1 en valeur linéaire) ce qui permet d'éviter le clipping. Cela est nécessaire car c'est l'amplitude dans le domaine fréquentiel qui compte et pas celle du domaine temporel (normaliser le pic absolu du signal dans le domaine temporel ne sert strictement à rien et peut être dangereux pour le matériel !). Le signal renvoyé est le signal normalisé dans le domaine temporel.

• readFile.m => function result = readFile(fileName)

Fonction qui permet de lire une impulsion au format .txt ascii (données brutes, les commentaires ne sont pas supportés !)

• simple_resize => function result=simple_resize(signal, newLength)

Fonction qui permet de redimensionner un signal, en l'agrandissant uniquement. Prend un signal et une longueur cible. Si la longueur cible est supérieure à la longueur actuelle, alors la fonction renvoie le signal avec un silence ajouté à la fin (des valeurs 0) de sorte que le signal retourné soit de la nouvelle taille souhaitée.

• trunc_size.m => function result=trunc_size(signal, newLength, window)

Fonction qui prend un signal dans le domaine temporel, une longueur cible, et (optionel) une fenêtre. Le signal est tronqué autour de son pic maximal en valeur absolue, puis fenêtré par la fenêtre si spécifié. La seule fenêtre disponible actuellement est la "blackman", c'est celle que j'utilise pour toutes mes troncatures.

Fonctionnement / utilisation

Comme vous l'aurez compris, la partie principale du programme se trouve dans le fichier "generate_one.m" qui contient le code pour lancer la génération d'une correction, et "generate_correction.m", qui correspond au code pour le calcul de cette correction.

Le contenu principal de "generate_one.m" est le suivant :

measure=readFile("data/sub.txt");
target=readFile("data/target_sub.txt");

Ces 2 lignes appellent la fonction readFile, qui lit des données brutes contenues dans un fichier txt. Ici il s'agit de charger dans la variable "measure" et "target" les impulsions de mesure et cible nécessaires au calcul de la correction. Vous êtes libre d'adapter ces lignes à vos besoins ! Par exemple si vous souhaitez lire un fichier .wav contenant l'impulsion, vous pourriez utiliser la ligne suivante :

measure=wavread("data/sub.wav")

Ce qui aurait pour effet de lire le fichier .wav et de le charger dans la même variable. Ce qui importe, c'est que la variable "measure" contienne votre impulsion.

Ensuite la ligne suivante appelle la fonction "generate_correction" en lui passant en paramètres la mesure et la courbe cible, et la correction obtenue et enregistrée dans la variable "correction" :

correction=generate_correction(measure, target);

Enfin il s'agit d'exporter les données que contient cette variable sous forme de fichier .txt ou .wav pour pouvoir faire la convolution avec le logiciel de son choix (pour ma part j'utilise le .txt pour BruteFir et le .wav pour ConvolverVST) :

En fichier .txt brut (sans commentaires) :

save "output/correction_sub.txt" correction -ascii;

En fichier .wav 32 bits a 48khz :

wavwrite(correction, 48000, 32, "output/sub.txt");

Voilà rien de bien compliqué jusque là ! Pour lancer la génération de la correction il suffit donc, depuis Octave (et dans le répertoire du logiciel), de taper :

generate_one;

et de valider par [Entrée].

Le code contenu dans le fichier generate_one.m sera donc lancé, chargera les impulsions, générera la correction, et l'exportera en .txt ou .wav !

Enfin vous pouvez éditer le fichier generate_correction.m pour voir comment est effectué la correction et quels sont les paramètres que j'utilise. J'ai en premier lieu écrit un bloc de lignes permettant de "régler" les étapes de la correction. Je conseille de laisser toutes les étapes activées (sur 1), vous pouvez les désactiver pour comprendre le fonctionnement de la correction par exemple, en changeant les valeurs 1 par 0.

constantSignalTrick=1 # This avoid low frequency error when calculating target/measure

Lorsque j'appliquais l'inversion courbe cible / measure, j'obtenais une correction contenant des basses fréquences d'une amplitude non négligeable (jusqu'à -10dB !). J'ai tenté de rajouter du "bruit" au signal, en vain car l'énergie du bruit est répartie avant/après l'impulsion (étalée équitablement sur tout le domaine temporel) et la correction générée, bien que parfaite dans le domaine fréquentiel en amplitude/phase, n'est pas du tout bonne en terme temporel (visible particulièrement sur le spectrogramme ou sur la forme d'impulsion dans le domaine temporel et en échelle logarithmique). La meilleure solution que j'ai trouvé et qui semble fonctionner parfaitement est de rajouter une composante continue au signal de mesure. Alors l'inversion ne comporte plus de basse fréquence et la réponse est bonne dans le domaine fréquenciel et temporel cette fois.

truncatingMeasure=1 # Because it's useless in general to have a very long impulse, I truncate it to a

Cette étape permet de tronquer et fenêtrer la mesure à +/-10000 samples (longueur de 20000 samples). Vous pouvez désactiver cette étape ou changer la longueur de troncature plus loin dans le fichier. Il est en effet inutile de travailler sur un signal trop long, et le fait de tronquer la mesure permet de diminuer l'énergie inutile dû au bruit. En tronquant on augmente donc la rapidité de calcul mais surtout le rapport SNR donc on diminue la THD.

correctionCleaningFiltering=1

Comme son nom l'indique cette étape permet de faire un "nettoyage" de la correction obtenue. En effet en raison du bruit (bruit numérique et bruit de mesure), le signal de correction obtenue présente toujours de l'énergie inutile dans le bas du spectre et, plus particulièrement, le haut du spectre. Pour détecter la fréquence de nettoyage passe-haut/passe-bas, je me base sur la courbe cible, et je descend dans le spectre jusqu'à obtenir une amplitude inférieure à une certaine dynamique. Par défaut cette valeur de dynamique est fixée à -80dB, vous pouvez la changer dans le fichier un peu plus loin. Lorsque la dynamique vaut -80dB, je considère qu'à partir de cette fréquence il faut supprimer tout ce qui se trouve en dessous. Idem pour le haut du spectre. Je génère alors un filtre passe-bande brickwall (pente maximale) pour supprimer tout ce qui se trouve dans dehors de ces 2 fréquences.

correctionTruncatingWindowing=1

La correction obtenue est finalement tronquée et fenêtrée par un fenêtrage blackman.

correctionNormalizeFreqPeak=1

Enfin je normalise le signal pour que le pic dans le domaine fréquenciel ne dépasse pas 1 (0dBFS) pour éviter tout clipping possible.

Pour mon besoin personnel, effectuant une correction sur 2x3 voies + subwoofer + 2 large bande arrière, j'utilise un fichier "generate_all.m" pour générer tout les corrections en un seul script. Les fichiers d'impulsions sont dans le dossier "data" et les corrections sont générées dans le dossier "output". Il me reste plus qu'à récupérer les corrections et les charger dans mon logiciel de convolution pour vérifier le résultat.

J'effectue alors une mesure en champ proche pour vérifier que la correction fonctionne comme souhaitée (réponse en fréquence en amplitude/phase correcte et SURTOUT spectrogramme sans résonance anormale ni énergie avant le pic d'impulsion).

Enfin j'effectue le calage des voies au niveau temporel en ajustant les délais pour chacune d'elles, puis je corrige la réponse fréquencielle en amplitude avec un égaliseur IIR pour compenser les erreurs diverses apportées par les résonances de pièce et autres. A partir de là, la correction est terminée !