Ich möchte eine Rauschformung in einer 100-kHz-16-Bit-Anwendung durchführen, um das gesamte Quantisierungsrauschen auf das 25-kHz-50-kHz-Band mit minimalem Rauschen im DC-25-kHz-Band zu verschieben.

Ich habe mathematica so eingerichtet, dass durch Verstärkungslernen ein Fehlerfilterkern mit 31 Stichproben erstellt wird, der gut funktioniert: Nach ein wenig Lernen kann ich das Hochfrequenzrauschen um ca. 16 dB verbessern, wobei das Niederfrequenzband (das Mittellinie ist der ungeformte Dither-Geräuschpegel. Dies steht im Einklang mit dem Rauschformungssatz "Gerzon-Craven".

Nun zu meinem Problem:

Ich kann es nicht schaffen, das Rauschen selbst nach ausgiebigem Lernen noch mehr zu formen, obwohl das Gerzon-Craven-Theorem es nicht verbietet. Zum Beispiel sollte es möglich sein, eine Reduzierung von 40 dB im niedrigen Band und eine Verbesserung von 40 dB im hohen Band zu erreichen.

Gibt es also eine andere grundlegende Grenze, auf die ich stoße?

Ich habe versucht, die Shannon-Rausch- / Abtast- / Informationssätze zu untersuchen, aber nachdem ich eine Weile damit herumgespielt hatte, gelang es mir nur, eine einzige Grenze daraus abzuleiten: den Gerzon-Craven-Satz, der ein direktes Ergebnis des Shannon-Satzes zu sein scheint.

Jede Hilfe wird geschätzt.

EDIT: mehr Infos

Zunächst der Filterkern, der die obige Rauschformung erzeugt. Beachten Sie, dass sich das neueste Sample auf der rechten Seite befindet. Die numerischen Werte des BarChart sind auf 0,01 gerundet: {-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} (Nicht genau der Balken, erzeugt jedoch eine ähnliche Kurve )

Ein weiterer Hinweis zur Implementierung von Fehlerrückmeldungen:

Ich habe zwei verschiedene Implementierungen des Fehlerfeedbacks ausprobiert. Zuerst habe ich die gerundete Ausgangsstichprobe mit dem gewünschten Wert verglichen und diese Abweichung als Fehler verwendet. Zweitens habe ich die gerundete Ausgangsstichprobe mit der (Eingabe + Fehlerrückmeldung) verglichen. Obwohl beide Methoden ganz unterschiedliche Kernel erzeugen, scheinen sich beide mit etwa der gleichen Intensität der Rauschformung abzugleichen. Die hier veröffentlichten Daten verwenden die zweite Implementierung.

Hier ist der Code, der zur Berechnung der digitalisierten Wellenabtastwerte verwendet wird. Schritt ist die Schrittgröße zum Runden. Welle ist die nicht digitalisierte Wellenform (normalerweise nur Nullen, wenn kein Signal angelegt wird).

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
  ]

Die Bewehrungsmethode:

Die "Punktzahl" wird berechnet, indem das Rauschleistungsspektrum betrachtet wird. Ziel ist es, die Rauschleistung im Band DC-25kHz zu minimieren. Ich bestrafe kein Rauschen im Hochfrequenzband, so dass ein willkürlich hohes Rauschen die Punktzahl nicht verringern würde. Ich führe Rauschen in die Kernelgewichte ein, um zu lernen. Vielleicht bin ich deshalb in einem (sehr breiten und tiefen) lokalen Minimum, aber ich halte dies für äußerst unwahrscheinlich.

Vergleich zum Standardfilterdesign:

Mit Mathematica können Filter iterativ generiert werden. Diese können einen viel besseren Kontrast als 36 dB haben, wenn ihr Frequenzgang aufgezeichnet wird; bis zu 80-100 dB. Numerische Werte: {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}

Wenn jedoch diejenigen in der tatsächlichen Rauschformung angewendet werden, werden sie (a) auf den gleichen Kontrast von ~ 40 dB geklemmt, (b) schlechter als das gelernte Filter, das tatsächlich keine Rauschdämpfung ausführt.

Grundlegendes Dithering ohne Noise Shaping

Die grundlegende Dither-Quantisierung ohne Rauschformung funktioniert folgendermaßen:

Abbildung 1. Grundlegendes Diagramm des Dither-Quantisierungssystems. Rauschen ist ein dreieckiges Dithering mit dem Mittelwert Null mit einem maximalen Absolutwert von 1. Die Rundung erfolgt auf die nächste Ganzzahl. Der Restfehler ist die Differenz zwischen Ausgabe und Eingabe und wird nur für die Analyse berechnet.

Das dreieckige Zittern erhöht die Varianz des resultierenden Restfehlers um den Faktor 3 (von bis ), entkoppelt jedoch den Mittelwert und die Varianz des Nettoquantisierungsfehlers vom Wert des Eingangssignals. Das bedeutet, dass das Nettofehlersignal nicht mit dem Eingang korreliert ist, aber höhere Momente nicht entkoppelt sind, so dass es kein wirklich völlig unabhängiger Zufallsfehler ist, aber niemand hat festgestellt, dass Menschen eine Abhängigkeit der höheren Momente im Nettofehlersignal von der hören können Eingangssignal in einer Audioanwendung.$\frac{1}{12}$$\frac{1}{4}$

Mit einem unabhängigen additiven Restfehler hätten wir ein einfacheres Modell des Systems:

Abbildung 2. Approximation der grundlegenden Dither-Quantisierung. Restfehler ist weißes Rauschen.

Im ungefähren Modell wird die Ausgabe einfach plus unabhängiger Restfehler des weißen Rauschens eingegeben.

Dithering mit Noise Shaping

Ich kann Mathematica nicht sehr gut lesen, daher werde ich anstelle Ihres Systems das System von Lipshitz et al. " Minimal hörbare Geräuschformung " J. Audio Eng. Soc., Band 39, Nr. 11, November 1991:

Abbildung 3. Lipshitz et al. Systemdiagramm von 1991 (nach Abb. 1). Der Filter (im Text kursiv geschrieben) enthält eine Verzögerung von einem Abtastwert, sodass er als Fehlerrückkopplungsfilter verwendet werden kann. Lärm ist dreieckiges Zittern.

Wenn der Restfehler unabhängig von aktuellen und vergangenen Werten von Signal A ist, haben wir ein einfacheres System:

Figure 4. Ein ungefähres Modell von Lipshitz et al. 1991 System. Der Filter ist der gleiche wie in Fig. 3 und enthält eine Verzögerung von einem Abtastwert. Es wird nicht mehr als Rückkopplungsfilter verwendet. Restfehler ist weißes Rauschen.

In dieser Antwort werde ich mit dem leichter zu analysierenden Näherungsmodell arbeiten (Abb. 4). Im Original haben Lipshitz et al. 1991 System hat Filter eine generische Filterform mit unendlicher Impulsantwort (IIR), die sowohl IIR- als auch FIR-Filter (Finite Impulse Response) abdeckt. Im Folgenden nehmen wir an, dass Filter ein FIR-Filter ist, da ich aufgrund meiner Experimente mit Ihren Koeffizienten glaube, dass dies das ist, was Sie in Ihrem System haben. Die Übertragungsfunktion von Filter ist:

Der Faktor repräsentiert eine Verzögerung von einer Abtastung. Im Näherungsmodell gibt es auch einen direkten Summierungspfad, der vom Restfehler ausgegeben wird. Dies wird mit der negierten Ausgabe von Filter summiert und bildet die vollständige Rauschformungsfilter-Übertragungsfunktion:$z^{-1}$

Um von Ihren Filterkoeffizienten , die Sie in der Reihenfolge , zu den Polynomkoeffizienten , dem Vorzeichen der Koeffizienten wird geändert, um die Negation der Filterausgabe im Systemdiagramm zu berücksichtigen , und der Koeffizient wird an das Ende angehängt (bis im Octave-Skript unten), und schließlich wird die Liste umgekehrt (bis ): 1 , b 1 , b 2 , b 3 , … b 0 = $\ldots, -b_3, -b_2, -b_1$$1, b_1, b_2, b_3, \ldots$$b_0 = 1$horzcatflip

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)

Das Skript zeichnet den Größenfrequenzgang und die Nullstellen des vollständigen Rauschformungsfilters auf:

Abbildung 5. Größenfrequenzgang des vollständigen Rauschformungsfilters.

Abbildung 6. Z-Ebenen-Diagramm der Pole ( ) und Nullen ( ) des Filters. Alle Nullen befinden sich innerhalb des Einheitskreises, sodass das vollständige Rauschformungsfilter eine minimale Phase aufweist.$\times$$\circ$

Ich denke, das Problem des Findens der Filterkoeffizienten kann als das Problem des Entwurfs eines Minimalphasenfilters mit einem führenden Koeffizienten von 1 umformuliert werden. Wenn der Frequenzgang solcher Filter inhärente Einschränkungen aufweist, werden diese Einschränkungen auf äquivalente Einschränkungen übertragen bei der Rauschformung, die solche Filter verwendet.

Umstellung vom Allpoldesign auf Minimum-Phase-FIR

Ein Verfahren zum Entwurf verschiedener, aber in vielerlei Hinsicht äquivalenter Filter ist in Stojanović et al. , "Allpoliges rekursives digitales Filterdesign basierend auf Ultraschallpolynomen", Radioengineering, Band 23, Nr. 3, September 2014. Sie berechnen Nennerkoeffizienten der Übertragungsfunktion eines allpoligen IIR-Tiefpassfilters. Diese haben immer einen führenden Nennerkoeffizienten von 1 und alle Pole innerhalb des Einheitskreises, was stabile IIR-Filter erfordert. Wenn diese Koeffizienten als Koeffizienten des FIR-Rauschformungsfilters mit minimaler Phase verwendet werden, ergeben sie im Vergleich zum Tiefpass-IIR-Filter einen invertierten Hochpass-Frequenzgang (Übertragungsfunktions-Nennerkoeffizienten werden zu Zählerkoeffizienten). In Ihrer Notation befindet sich ein Satz von Koeffizienten aus diesem Artikel {-0.0076120, 0.0960380, -0.5454670, 1.8298040, -3.9884220, 5.8308660, -5.6495140, 3.3816780}, der für die Rauschformungsanwendung getestet werden könnte, obwohl er nicht genau der Spezifikation entspricht:

Figure 7. Größenfrequenzgang des FIR-Filters unter Verwendung von Koeffizienten von Stojanović et al. 2014.

Figure 8. Pol-Null-Diagramm des FIR-Filters unter Verwendung von Koeffizienten von Stojanović et al. 2014.

Die allpolige Übertragungsfunktion ist:

Sie können also ein stabiles allpoliges IIR-Tiefpassfilter entwerfen und die Koeffizienten als Koeffizienten verwenden, um ein Hochphasen-FIR-Filter mit minimaler Phase und einem führenden Koeffizienten von 1 zu erhalten.$a$$b$

Um ein Allpolfilter zu entwerfen und dieses in ein FIR-Filter mit minimaler Phase umzuwandeln, können Sie keine IIR-Filterentwurfsmethoden verwenden, die von einem analogen Prototypfilter ausgehen und die Pole und Nullen mithilfe der bilinearen Transformation in den digitalen Bereich abbilden . Dazu gehören cheby1, cheby2und ellipin Octave und Pythons SciPy. Diese Methoden geben Nullen vom Ursprung der Z-Ebene weg, sodass der Filter nicht vom erforderlichen Allpoltyp ist.

Antwort auf die theoretische Frage

Wenn es Ihnen egal ist, wie viel Rauschen bei Frequenzen über einem Viertel der Abtastfrequenz auftritt, dann haben Lipshitz et al. 1991 geht Ihre Frage direkt an:

Für solche Gewichtungsfunktionen, die über einen Teil des Bandes auf Null gehen, gibt es keine theoretische Grenze für die gewichtete Rauschleistungsreduzierung, die aus der Schaltung von Fig. 1 erhalten werden kann. Dies wäre der Fall, wenn man beispielsweise annimmt, dass die Das Ohr hat eine Empfindlichkeit von beispielsweise zwischen 20 kHz und der Nyquist-Frequenz von Null und wählt die Gewichtungsfunktion, um diese Tatsache widerzuspiegeln.

Ihre Abbildung 1 zeigt einen Rauschformer mit einer generischen IIR-Filterstruktur mit sowohl Polen als auch Nullen, die sich so stark von der FIR-Struktur unterscheidet, die Sie derzeit haben. Was sie sagen, gilt jedoch auch dafür, da eine FIR-Filter-Impulsantwort sein kann willkürlich nahe an der Impulsantwort eines gegebenen stabilen IIR-Filters gemacht.

Oktavskript für Filterdesign

Hier ist ein Oktavskript zur Koeffizientenberechnung mit einer anderen Methode, die meiner Meinung nach der von Stojanovici et al. 2014 Methode parametrisiert als mit der richtigen Wahl meines Parameters.$\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

Es beginnt mit einem Dolph-Chebyshev-Fenster, während die Koeffizienten es mit sich selbst falten, um die Nullen der Übertragungsfunktion zu verdoppeln, und fügt dem mittleren Abgriff eine Zahl hinzu, die den Frequenzgang "anhebt" (wobei der mittlere Abgriff als Nullzeit betrachtet wird) dass es überall positiv ist, die Nullen findet, Nullen entfernt, die außerhalb des Einheitskreises liegen, die Nullen wieder in Koeffizienten umwandelt (der führende Koeffizient von polyist immer 1) und das Vorzeichen jedes zweiten Koeffizienten umdreht, um den Filter hochpass zu machen . Die Ergebnisse des Skripts (eine ältere, aber fast gleichwertige Version) sehen vielversprechend aus:

Abbildung 9. Größenfrequenzgang des Filters aus (einer älteren, aber fast gleichwertigen Version) des obigen Skripts.

Abbildung 10. Pol-Null-Diagramm des Filters aus (einer älteren, aber fast gleichwertigen Version) des obigen Skripts.

Die Koeffizienten aus (einer älteren, aber fast äquivalenten Version von) dem obigen Skript in Ihrer 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}. Die Zahlen sind groß, was zu numerischen Problemen führen kann.

Oktavimplementierung der Rauschformung

Schließlich habe ich die Rauschformung in Octave selbst implementiert und bekomme keine Probleme wie Sie. Aufgrund unserer Diskussion in Kommentaren bestand die Einschränkung in Ihrer Implementierung meiner Meinung nach darin, dass das Rauschspektrum unter Verwendung eines rechteckigen Fensters, auch "kein Fenster" genannt, bewertet wurde , das das Hochfrequenzspektrum auf die niedrigen Frequenzen verschüttete.

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. Quantisierungsrauschspektralanalyse aus der obigen Oktavimplementierung der Rauschformung für ein Eingangssignal mit konstantem Nullpunkt. Horizontale Achse: Normalisierte Frequenz. Schwarz: keine Rauschformung ( c = [1];), Rot: Ihr Originalfilter, Blau: Der Filter aus Abschnitt "Oktavskript für Filterdesign".

Abbildung 12. Zeitbereichsausgabe aus der obigen Oktavimplementierung der Rauschformung für ein Eingangssignal mit konstantem Nullpunkt. Horizontale Achse: Probennummer, vertikale Achse: Probenwert. Rot: Ihr Originalfilter, Blau: Der Filter aus Abschnitt "Oktavskript für Filterdesign".

Das extremere Rauschformungsfilter (blau) führt zu sehr großen quantisierten Ausgangsabtastwerten, selbst bei Null-Eingang.