Je veux faire une mise en forme du bruit dans une application 100 kHz, 16 bits, de manière à déplacer tout le bruit de quantification vers la bande 25 kHz-50 kHz, avec un bruit minimal dans la bande DC-25 kHz.

J'ai mis en place mathématique pour créer un noyau de filtre d'erreur à 31 échantillons via un apprentissage par renforcement qui fonctionne bien: après un peu d'apprentissage, je peux obtenir une amélioration d'environ 16 dB du bruit haute fréquence pour environ la même quantité de réduction dans la bande basse fréquence (la la ligne centrale est le niveau de bruit de vibration non formé). Ceci est conforme au théorème de mise en forme du bruit "Gerzon-Craven".

Passons maintenant à mon problème:

Je n'arrive pas à façonner le bruit encore plus même après un apprentissage approfondi, bien que le théorème de Gerzon-Craven ne l'interdise pas. Par exemple, il devrait être possible d'obtenir une réduction de 40 dB dans la bande basse et une amélioration de 40 dB dans la bande haute.

Y a-t-il donc une autre limite fondamentale que je rencontre?

J'ai essayé d'examiner les théorèmes de bruit / échantillonnage / information de Shannon, mais après y avoir longuement réfléchi, je n'ai réussi à en déduire qu'une seule limite: le théorème de Gerzon-Craven, qui semble être un résultat direct du théorème de Shannon.

Toute aide est appréciée.

EDIT: plus d'infos

Tout d'abord, le noyau de filtre qui produit la mise en forme du bruit ci-dessus, notez que l' échantillon le plus récent est sur le côté droit . Valeurs numériques du graphique à barres arrondies à 0,01: {-0,16, 0,51, -0,74, 0,52, -0,04, -0,25, 0,22, -0,11, -0,02, 0,31, -0,56, 0,45, -0,13, 0,04, -0,14, 0,12, -0,06, 0,19, -0,22, -0,15, 0,4, 0,01, -0,41, -0,1, 0,84, -0,42, -0,81, 0,91, 0,75, -2,37, 2,29} (pas exactement le bar bar mais produit une courbe similaire )

Une autre note sur la mise en œuvre du retour d'erreur:

J'ai essayé deux implémentations différentes du retour d'erreur. J'ai d'abord comparé l'échantillon de sortie arrondi à la valeur souhaitée et j'ai utilisé cet écart comme erreur. Deuxièmement, j'ai comparé l'échantillon de sortie arrondi à (entrée + retour d'erreur). Bien que les deux méthodes produisent des noyaux assez différents, les deux semblent se stabiliser à environ la même intensité de mise en forme du bruit. Les données publiées ici utilisent la deuxième implémentation.

Voici le code utilisé pour calculer les échantillons d'ondes numérisés. étape est la taille d'étape pour l'arrondi. wave est la forme d'onde non numérisée (généralement des zéros uniquement quand aucun signal n'est appliqué).

TestWave[kernel_?VectorQ] := 
 Module[{k = kernel, nf, dith, signals, twave, deltas},
  nf = Length@k;
  dith = RandomVariate[TriangularDistribution[{-1, 1}*step], l];
  signals = deltas = Table[0, {l}];
  twave = wave;
  Do[
   twave[[i]] -= k.PadLeft[deltas[[;; i - 1]], nf];
   signals[[i]] = Round[twave[[i]] + dith[[i]], step];
   deltas[[i]] = signals[[i]] - twave[[i]];
   , {i, l}];
  signals
  ]

La méthode de renforcement:

Le "score" est calculé en examinant le spectre de puissance de bruit. L'objectif est de minimiser la puissance de bruit dans la bande DC-25kHz. Je ne pénalise pas le bruit dans la bande des hautes fréquences, donc un bruit arbitrairement élevé ne diminuerait pas le score. J'introduis du bruit dans les poids du noyau pour apprendre. Peut-être, par conséquent, je suis dans un minimum local (très large et profond), mais je considère cela extrêmement peu probable.

Comparaison avec la conception de filtre standard:

Mathematica permet de générer des filtres de manière itérative. Ceux-ci peuvent avoir un contraste bien meilleur que 36 dB lorsque leur réponse en fréquence est tracée; jusqu'à 80-100 dB. Valeurs numériques: {0,024, -0,061, -0,048, 0,38, -0,36, -0,808, 2,09, -0,331, -4,796, 6,142, 3,918, -17,773, 11,245, 30,613, -87,072, 113,676, -87,072, 30,613, 11,245 , -17,773, 3,918, 6,142, -4,796, -0,331, 2,09, -0,808, -0,36, 0,38, -0,048, -0,061, 0,024}

Cependant, lors de l'application de ceux dans la mise en forme réelle du bruit, ils sont (a) fixés au même contraste de ~ 40 dB, (b) fonctionnent moins bien que le filtre appris sans aucune atténuation du bruit.

noise tobalt
la source

+1, question très intéressante. Avez-vous essayé d'augmenter l'ordre du filtre au-dessus de 31 pressions? La suppression de 40 dB semble un peu élevée pour un FIR à 31 prises.

A_A

@Olli, je ne crois pas comprendre complètement. Je peux poster le noyau du filtre si c'est ce qui vous intéresse. En termes simples, il existe des poids oscillatoires qui forcent l'erreur à l'alternative -> la déplace vers les hautes fréquences.

2018 à 10h08

@tobalt de la conception de filtre "classique", il est attendu que les filtres plus longs soient plus raides et / ou aient plus d'atténuation dans la bande d'arrêt et / ou aient moins d'ondulation dans la bande passante. Maintenant, je suppose que votre méthode de renforcement récompense la pente plus que l'atténuation après un certain point; quelle est la méthode que vous utilisez pour renforcer?

Marcus Müller

Vous voudrez peut-être consulter la section Conception de filtre de Mathematica. Vous pouvez peut-être simplement définir les spécifications de votre filtre et utiliser l'une des techniques existantes pour renvoyer un filtre qui les satisfait.

A_A

C'est donc définitivement (éventuellement itératif) la conception du filtre. Obtenez vos spécifications de filtre (exactement comme vous les avez publiées ici) et essayez de créer un filtre via cette fonction (la plus simple de son type) et voyez ce qu'il propose. Ce serait bien de voir les coefficients de cette fonction par rapport à ceux avec lesquels l'apprentissage de renforcement revient. Notez également le type d'ordre du filtre qu'il propose, je suppose qu'il serait supérieur à 31. Doit-il être "adaptatif" au signal en passant?

A_A

Tramage de base sans mise en forme du bruit

La quantification de base tramée sans mise en forme du bruit fonctionne comme ceci:

Figure 1. Schéma du système de quantification de base tramé. Le bruit est un tramage triangulaire à moyenne nulle avec une valeur absolue maximale de 1. L'arrondi est à l'entier le plus proche. L'erreur résiduelle est la différence entre la sortie et l'entrée, et est calculée pour l'analyse uniquement.

$\frac{1}{12}$ $\frac{1}{4}$

Avec une erreur résiduelle additive indépendante, nous aurions un modèle plus simple du système:

Figure 2. Approximation de la quantification de base tramée. L'erreur résiduelle est le bruit blanc.

Dans le modèle approximatif, la sortie est simplement entrée plus une erreur résiduelle de bruit blanc indépendante.

Tramage avec mise en forme du bruit

Je ne peux pas très bien lire Mathematica, donc au lieu de votre système, j'analyserai le système de Lipshitz et al. " Mise en forme du bruit minimalement audible " J. Audio Eng. Soc., Vol.39, No.11, novembre 1991:

Figure 3. Lipshitz et al. Schéma du système de 1991 (adapté de leur Fig. 1). Le filtre (en italique dans le texte) comprend un délai d'un échantillon afin qu'il puisse être utilisé comme filtre de retour d'erreur. Le bruit est un tramage triangulaire.

Si l'erreur résiduelle est indépendante des valeurs actuelles et passées du signal A, nous avons un système plus simple:

Figure 4. Un modèle approximatif de Lipshitz et al. Système de 1991. Le filtre est le même que sur la figure 3 et comprend un retard d'un échantillon. Il n'est plus utilisé comme filtre de rétroaction. L'erreur résiduelle est le bruit blanc.

Dans cette réponse, je travaillerai avec le modèle approximatif plus facilement analysable (Fig. 4). Dans l'original Lipshitz et al. Système de 1991, Filter a une forme de filtre générique à réponse impulsionnelle infinie (IIR) qui couvre à la fois les filtres IIR et à réponse impulsionnelle finie (FIR). Dans ce qui suit, nous supposerons que Filter est un filtre FIR, comme je le crois d'après mes expériences avec vos coefficients, c'est ce que vous avez dans votre système. La fonction de transfert du filtre est:

H_{F i l t e r} (z) = - b_{1} z^{- 1} - b_{2} z^{- 2} - b_{3} z^{- 3} - \dots

$H_\mathit{Filter}(z) = -b_1z^{-1} - b_2z^{-2} - b_3z^{-3} - \ldots$

$z^{-1}$

H (z) = 1 - H_{F i l t e r} (z) = 1 + b_{1} z^{- 1} + b_{2} z^{- 2} + b_{3} z^{- 3} + \dots .

$H(z) = 1 - H_\mathit{Filter}(z) = 1 + b_1z^{-1} + b_2z^{-2} + b_3z^{-3} + \ldots.$

$\ldots, -b_3, -b_2, -b_1$ $1, b_1, b_2, b_3, \ldots$ $b_0 = 1$ horzcatdans le script Octave ci-dessous), et enfin la liste est inversée (par flip):

pkg load signal
b = [-0.16, 0.51, -0.74, 0.52, -0.04, -0.25, 0.22, -0.11, -0.02, 0.31, -0.56, 0.45, -0.13, 0.04, -0.14, 0.12, -0.06, 0.19, -0.22, -0.15, 0.4, 0.01, -0.41, -0.1, 0.84, -0.42, -0.81, 0.91, 0.75, -2.37, 2.29];
c = flip(horzcat(-b, 1));
freqz(c)
zplane(c)

Le script trace la réponse en fréquence d'amplitude et les emplacements zéro du filtre de mise en forme du bruit complet:

Figure 5. Réponse en fréquence d'amplitude du filtre de mise en forme du bruit complet.

$\times$ $\circ$

Je pense que le problème de trouver les coefficients du filtre peut être reformulé comme le problème de la conception d'un filtre à phase minimale avec un coefficient de tête de 1. S'il existe des limitations inhérentes à la réponse en fréquence de ces filtres, alors ces limitations sont transférées à des limitations équivalentes dans la mise en forme du bruit qui utilise de tels filtres.

Conversion de la conception omnipolaire en FIR à phase minimale

Une procédure de conception de filtres différents mais à bien des égards équivalents est décrite dans Stojanović et al. , "Conception de filtres numériques récursifs tous pôles basée sur des polynômes ultraspheriques", Radioengineering, vol 23, no 3, septembre 2014. Ils calculent les coefficients de dénominateur de la fonction de transfert d'un filtre passe-bas multipolaire IIR. Ceux-ci ont toujours un coefficient de dénominateur avancé de 1 et ont tous les pôles à l'intérieur du cercle unitaire, une exigence de filtres IIR stables. Si ces coefficients sont utilisés comme coefficients du filtre de mise en forme du bruit FIR à phase minimale, ils donneront une réponse en fréquence passe-haut inversée par rapport au filtre passe-bas IIR (les coefficients du dénominateur de la fonction de transfert devenant des coefficients du numérateur). Dans votre notation, un ensemble de coefficients de cet article est {-0.0076120, 0.0960380, -0.5454670, 1.8298040, -3.9884220, 5.8308660, -5.6495140, 3.3816780}, qui pourrait être testé pour l'application de mise en forme du bruit, bien qu'il ne soit pas exactement conforme aux spécifications:

Figure 7. Réponse en fréquence d'amplitude du filtre FIR à l'aide des coefficients de Stojanović et al. 2014.

Figure 8. Diagramme à pôles zéro du filtre FIR à l'aide des coefficients de Stojanović et al. 2014.

La fonction de transfert omnipolaire est:

H (z) = \frac{1}{1 + a_{1} z^{- 1} + a_{2} z^{- 2} + a_{3} z^{- 3} + \dots}

$H(z) = \frac{1}{1 + a_1z^{-1} + a_2z^{-2} + a_3z^{-3} + \ldots}$

$a$ $b$

Pour concevoir un filtre omnipolaire et le convertir en un filtre FIR à phase minimale, vous ne pourrez pas utiliser les méthodes de conception de filtre IIR qui partent d'un filtre prototype analogique et mappent les pôles et les zéros dans le domaine numérique à l'aide de la transformation bilinéaire. . Cela inclut cheby1, cheby2et ellipdans Octave et Scyth de Python. Ces méthodes éloigneront les zéros de l'origine du plan z, de sorte que le filtre ne sera pas du type multipolaire requis.

Réponse à la question théorique

Si vous ne vous souciez pas de la quantité de bruit à des fréquences supérieures au quart de la fréquence d'échantillonnage, alors Lipshitz et al. 1991 répond directement à votre question:

Pour de telles fonctions de pondération, qui vont à zéro sur une partie de la bande, il n'y a pas de limite théorique à la réduction de puissance de bruit pondérée pouvant être obtenue à partir du circuit de la figure 1. Ce serait le cas si, par exemple, on suppose que le l'oreille a une sensibilité nulle entre, disons, 20 kHz et la fréquence de Nyquist, et choisit la fonction de pondération pour refléter ce fait.

Leur Fig 1. montre un shaper de bruit avec une structure de filtre IIR générique avec des pôles et des zéros, si différente de la structure FIR que vous avez en ce moment, mais ce qu'ils disent s'applique également à cela, car une réponse impulsionnelle du filtre FIR peut être arbitrairement proche de la réponse impulsionnelle de tout filtre IIR stable donné.

Script d'octave pour la conception de filtres

$\nu=0$ dip

pkg load signal
N = 14; #number of taps including leading tap with coefficient 1
att = 97.5; #dB attenuation of Dolph-Chebyshev window, must be positive
dip = 2; #spectrum lift-up multiplier, must be above 1
c = chebwin(N, att);
c = conv(c, c);
c /= sum(c);
c(N) += dip*10^(-att/10);
r = roots(c);
j = (abs(r(:)) <= 1);
r = r(j);
c = real(poly(r));
c .*= (-1).^(0:(N-1)); #if this complains, then root finding has probably failed
freqz(c)
zplane(c)
printf('%f, ', flip(-c(2:end))), printf('\n'); #tobalt's format

Il commence par une fenêtre Dolph-Chebyshev comme coefficients, la convolution avec elle-même pour doubler les zéros de la fonction de transfert, ajoute au tapotement central un nombre qui "élève" la réponse en fréquence (considérant le tapotement médian comme étant à zéro) donc qu'il est partout positif, trouve les zéros, supprime les zéros qui sont en dehors du cercle unitaire, reconvertit les zéros en coefficients (le coefficient de tête de polyest toujours 1) et retourne le signe de chaque deuxième coefficient pour faire passer le filtre passe-haut . Les résultats (d'une version plus ancienne mais presque équivalente) du script semblent prometteurs:

Figure 9. Réponse en fréquence d'amplitude du filtre à partir d'une version plus ancienne mais presque équivalente du script ci-dessus.

Figure 10. Graphique à pôles zéro du filtre (d'une version plus ancienne mais presque équivalente) du script ci-dessus.

Les coefficients de (une ancienne mais presque version équivalente) , le script ci - dessus dans votre notation: {0.357662, -2.588396, 9.931419, -26.205448, 52.450624, -83.531276, 108.508775, -116.272581, 102.875781, -74.473956, 43.140431, -19.131434, 5.923468}. Les nombres sont importants, ce qui pourrait entraîner des problèmes numériques.

Implémentation octave de la mise en forme du bruit

Enfin, j'ai fait ma propre implémentation de la mise en forme du bruit dans Octave et je n'ai pas de problèmes comme vous. Sur la base de notre discussion dans les commentaires, je pense que la limitation dans votre mise en œuvre était que le spectre de bruit a été évalué en utilisant une fenêtre rectangulaire aka "pas de fenêtrage", qui a renversé le spectre haute fréquence vers les basses fréquences.

pkg load signal
N = length(c);
M = 16384; #signal length
input = zeros(M, 1);#sin(0.01*(1:M))*127;
er = zeros(M, 1);
output = zeros(M, 1);
for i = 1:M
  A = input(i) + er(i);
  output(i) = round(A + rand() - rand());
  for j = 2:N
    if (i + j - 1 <= M)
      er(i + j - 1) += (output(i) - A)*c(j);
    endif
  endfor
endfor
pwelch(output, max(nuttallwin(1024), 0), 'semilogy');

Figure 11. Analyse spectrale du bruit de quantification à partir de l'implémentation d'octave ci-dessus de la mise en forme du bruit pour un signal d'entrée à zéro constant. Axe horizontal: fréquence normalisée. Noir: pas de formation de bruit ( c = [1];), rouge: votre filtre d'origine, bleu: le filtre de la section "Script d'octave pour la conception du filtre".

Figure 12. Sortie dans le domaine temporel de l'implémentation d'octave ci-dessus de la mise en forme du bruit pour un signal d'entrée à zéro constant. Axe horizontal: numéro d'échantillon, axe vertical: valeur d'échantillon. Rouge: votre filtre d'origine, bleu: le filtre de la section "Script d'octave pour la conception du filtre".

Le filtre de mise en forme de bruit plus extrême (bleu) produit des valeurs d'échantillonnage de sortie quantifiées très importantes même pour une entrée nulle.

Olli Niemitalo
la source

@MattL. J'ai pensé à tort au début que le tobalt avait un filtre omnipolaire. J'ai réécrit ma réponse quand j'ai réalisé qu'il s'agissait d'un filtre FIR avec le premier coefficient 1. De plus, Gerzon-Craven aurait dit que le filtre doit être à phase minimale pour être optimal, et les coefficients optimisés de Tobalt donnent également un filtre à phase minimale. Ces exigences sont équivalentes à ce que les coefficients des filtres tout-pôles IIR ont donc je suggère d'emprunter des méthodes de conception à partir de là. Un IIR standard serait également une option.

Olli Niemitalo

J'ai isolé l'erreur: mon implémentation produit la même forme d'onde (dans le temps) que la vôtre. Cependant, la fonction Abs [Fourier [vague]] semble se heurter à un certain débordement / débordement interne, car le spectre renvoyé semble différent (étage supérieur)

tobalt

@Olli Niemitalo Ok, il semble que la FFT en octave utilise éventuellement le fenêtrage automatique? Après avoir appliqué une fenêtre Hann à la forme d'onde, je peux obtenir des FFT "correctes". Je vais tester brièvement l'intégrité de cette approche et éventuellement continuer à apprendre et publier le résultat. Merci pour tous vos efforts. J'ai marqué votre message comme réponse.