Chcę kształtować szumy w aplikacji o częstotliwości 100 kHz, 16 bitów, aby przenieść cały szum kwantyzacji na pasmo 25 kHz-50 kHz, przy minimalnym szumie w paśmie DC-25 kHz.

Skonfigurowałem matematykę, aby utworzyć 31-próbkowe jądro z filtrem błędów poprzez uczenie się wzmacniające, które działa dobrze: Po krótkim nauczeniu mogę uzyskać około ~ 16dB ulepszenia szumu wysokiej częstotliwości przy mniej więcej takiej samej redukcji w paśmie niskich częstotliwości ( linia środkowa to nieokształtowany poziom szumu ditherowego). Jest to zgodne z twierdzeniem o kształtowaniu hałasu „Gerzon-Craven”.

Teraz mój problem:

Nie potrafię jeszcze bardziej kształtować hałasu nawet po dogłębnym nauczeniu się, chociaż twierdzenie Gerzona-Cravena tego nie zabrania. Na przykład, powinno być możliwe osiągnięcie zmniejszenia 40 dB w dolnym paśmie i zwiększenie 40 dB w górnym paśmie.

Czy jest jeszcze inny podstawowy limit, na który wpadam?

Próbowałem przyjrzeć się twierdzeniom Shannona o hałasie / próbkowaniu / informacjach, ale po jakimś czasie manipulowania nim udało mi się wyprowadzić z tego tylko jedną granicę: twierdzenie Gerzona-Cravena, które wydaje się być bezpośrednim wynikiem twierdzenia Shannona.

Każda pomoc jest mile widziana.

EDYCJA: więcej informacji

Po pierwsze, jądro filtra, które wytwarza powyższe kształtowanie szumu, zauważ, że najnowsza próbka znajduje się po prawej stronie. Wartości liczbowe BarChart zaokrąglone do 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} (Niezupełnie charakterystyka słupkowa, ale daje podobną krzywą )

Kolejna uwaga na temat implementacji sprzężenia zwrotnego błędu:

Próbowałem dwóch różnych implementacji sprzężenia zwrotnego błędu. Najpierw porównałem zaokrągloną próbkę wyjściową do pożądanej wartości i wykorzystałem to odchylenie jako błąd. Po drugie porównałem zaokrągloną próbkę wyjściową do (wejście + informacja zwrotna o błędzie). Chociaż obie metody wytwarzają całkiem różne jądra, oba wydają się wyrównywać przy mniej więcej takiej samej intensywności kształtowania hałasu. Dane zamieszczone tutaj wykorzystują drugą implementację.

Oto kod używany do obliczania próbek fali cyfrowej. krok to wielkość kroku dla zaokrąglania. fala jest falą niezigitalizowaną (zazwyczaj tylko zera, gdy nie jest stosowany żaden sygnał).

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
  ]

Metoda zbrojenia:

„Wynik” jest obliczany na podstawie widma mocy szumu. Celem jest zminimalizowanie mocy szumu w paśmie DC-25kHz. Ja nie karania szum w wysokim paśmie częstotliwości, więc dowolnie wysoki hałas nie zmniejszyłyby wynik. Wprowadzam hałas do wag jądra do nauki. Może dlatego jestem w (bardzo szerokim i głębokim) minimum lokalnym, ale uważam to za niezwykle mało prawdopodobne.

Porównanie do standardowej konstrukcji filtra:

Mathematica pozwala na iteracyjne generowanie filtrów. Mogą mieć znacznie lepszy kontrast niż 36 dB, gdy wykreślone jest ich pasmo przenoszenia; do 80-100 dB. Wartości liczbowe: {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}

Jednakże, stosując te w rzeczywistym kształtowaniu szumu, są one (a) zaciśnięte na tym samym kontraście ~ 40dB, (b) działają gorzej niż wyuczony filtr faktycznie nie tłumiąc hałasu.

Podstawowe dithering bez kształtowania hałasu

Podstawowa kwantowa kwantyzacja bez kształtowania szumu działa w następujący sposób:

Ryc. 1. Podstawowy schemat systemu kwantyfikacji. Hałas jest zerową średnią trójkątną ditherem o maksymalnej wartości bezwzględnej 1. Zaokrąglenie jest do najbliższej liczby całkowitej. Błąd resztkowy jest różnicą między wyjściem a danymi wejściowymi i jest obliczany wyłącznie do analizy.

Trójkątny dither zwiększa wariancję wynikowego błędu resztkowego o współczynnik 3 (z do ), ale oddziela średnią i wariancję błędu kwantyzacji netto od wartości sygnału wejściowego. Oznacza to, że sygnał błędu netto nie jest skorelowany z wejściem, ale wyższe momenty nie są oddzielone, więc nie jest to tak naprawdę całkowicie niezależny błąd losowy, ale nikt nie ustalił, że ludzie mogą usłyszeć zależność wyższych momentów w sygnale błędu netto od sygnał wejściowy w aplikacji audio.$\frac{1}{12}$$\frac{1}{4}$

Przy niezależnym addytywnym błędzie resztkowym mielibyśmy prostszy model systemu:

Ryc. 2. Przybliżenie podstawowej kwantyfikacji skróconej. Błąd resztkowy to biały szum.

W przybliżonym modelu wyjście jest po prostu wprowadzane plus niezależny błąd resztkowy szumu białego.

Dithering z kształtowaniem hałasu

Nie potrafię dobrze czytać Mathematica, więc zamiast twojego systemu przeanalizuję system z Lipshitz i in. „ Minimalnie słyszalne kształtowanie hałasu ” J. Audio Eng. Soc., T. 39, nr 11, listopad 1991:

Ryc. 3. Lipshitz i in. Schemat systemu z 1991 r. (Na podstawie ich ryc. 1). Filtr (zaznaczony kursywą w tekście) zawiera w sobie jedno próbne opóźnienie, dzięki czemu można go wykorzystać jako filtr sprzężenia zwrotnego błędu. Hałas ma charakter trójkątny.

Jeśli błąd resztkowy jest niezależny od bieżących i przeszłych wartości sygnału A, mamy prostszy system:

Ryc. 4. Przybliżony model Lipshitza i in. System z 1991 r. Filtr jest taki sam jak na ryc. 3 i zawiera w sobie opóźnienie jednej próbki. Nie jest już używany jako filtr sprzężenia zwrotnego. Błąd resztkowy to biały szum.

W tej odpowiedzi będę pracować z łatwiejszym do analizy modelem przybliżonym (ryc. 4). W oryginalnej Lipshitz i in. System 1991, Filter ma ogólną formę filtra o nieskończonej odpowiedzi impulsowej (IIR), która obejmuje zarówno filtry IIR, jak i skończoną odpowiedź impulsową (FIR). W dalszej części założymy, że Filtr jest filtrem FIR, jak sądzę na podstawie moich eksperymentów z twoimi współczynnikami, że masz to w swoim systemie. Funkcja przesyłania filtra to:

Współczynnik reprezentuje opóźnienie jednej próby. W modelu przybliżonym istnieje również bezpośrednia ścieżka sumowania do wyjścia z błędu resztkowego. Jest to sumowane z zanegowanym wyjściem filtra , tworząc pełną funkcję przenoszenia filtra kształtującą szum:$z^{-1}$

Aby przejść od współczynników filtrowania , które podajesz w kolejności , do pełnej funkcji filtrowania kształtującego szumy współczynniki wielomianowe , znak współczynników zmienia się w celu uwzględnienia negacji danych wyjściowych filtru na schemacie systemowym, a współczynnik jest dołączany na końcu (przez w skrypcie Octave poniżej), a na koniec lista jest odwracana (o ):$\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)

Skrypt wykreśla pasmo przenoszenia wielkości i lokalizacje zerowe pełnego filtra kształtującego szum:

Rysunek 5. Pasmo przenoszenia wielkości filtra pełnego kształtującego szum.

Rysunek 6. Wykres Z w płaszczyźnie Z biegunów ( ) i zer ( ) filtra. Wszystkie zera znajdują się wewnątrz okręgu jednostki, więc pełny filtr kształtujący szum ma fazę minimalną.$\times$$\circ$

Myślę, że problem znalezienia współczynników filtra można przeformułować jako problem projektowania filtra fazy minimalnej o wiodącym współczynniku wynoszącym 1. Jeśli istnieją nieodłączne ograniczenia odpowiedzi częstotliwościowej takich filtrów, wówczas ograniczenia te przenoszone są na ograniczenia równoważne w kształtowaniu hałasu, który wykorzystuje takie filtry.

Konwersja z konstrukcji wielobiegunowej do FIR z fazą minimalną

Procedurę projektowania różnych, ale pod wieloma względami równoważnych filtrów opisano w Stojanović i in. , „Projektowanie wszystkich cyfrowych filtrów rekurencyjnych w oparciu o wielomiany ultradźwiękowe”, Radioengineering, tom 23, nr 3, wrzesień 2014. Obliczają współczynniki mianownika funkcji transferowej wielobiegunowego filtra dolnoprzepustowego IIR. Te zawsze mają wiodący współczynnik mianownika wynoszący 1 i mają wszystkie bieguny wewnątrz okręgu jednostki, co jest wymogiem stabilnych filtrów IIR. Jeżeli współczynniki te zostaną wykorzystane jako współczynniki filtra kształtującego szum FIR o minimalnej fazie, dadzą odwróconą górnoprzepustową odpowiedź częstotliwościową w porównaniu do dolnoprzepustowego filtra IIR (współczynniki mianownika funkcji przenoszenia stają się współczynnikami licznikowymi). W twojej notacji jest jeden zestaw współczynników z tego artykułu {-0.0076120, 0.0960380, -0.5454670, 1.8298040, -3.9884220, 5.8308660, -5.6495140, 3.3816780}, który można przetestować pod kątem aplikacji do kształtowania hałasu, chociaż nie jest on dokładnie zgodny ze specyfikacją:

Ryc. 7. Odpowiedź częstotliwościowa wielkości filtra FIR przy użyciu współczynników z Stojanović i in. 2014.

Rycina 8. Wykres bieguna zero dla filtra FIR przy użyciu współczynników z Stojanović i in. 2014.

Funkcja przenoszenia wszystkich biegunów to:

Tak więc, można zaprojektować stabilny wszystkich biegunów IIR filtr dolnoprzepustowy i użyć współczynników a współczynniki, aby uzyskać minimalną-fazowy filtr górnoprzepustowy FIR z wiodącym współczynnik 1.$a$$b$

Aby zaprojektować filtr wielobiegunowy i przekonwertować go na filtr FIR o minimalnej fazie, nie będzie można używać metod projektowania filtrów IIR, które zaczynają się od analogowego filtra prototypowego i odwzorowują bieguny i zera w domenie cyfrowej za pomocą transformacji dwuliniowej . Obejmuje to cheby1, cheby2oraz ellipw SciPy Octave'a i Pythona. Te metody podadzą zera od początku z płaszczyzny Z, więc filtr nie będzie wymaganego typu wielobiegunowego.

Odpowiedz na pytanie teoretyczne

Jeśli nie obchodzi cię, ile hałasu będzie na częstotliwościach przekraczających jedną czwartą częstotliwości próbkowania, to Lipshitz i in. 1991 odpowiada bezpośrednio na twoje pytanie:

Dla takich funkcji ważenia, które idą do zera na części pasma, nie ma teoretycznego ograniczenia ważonej redukcji mocy szumu możliwej do uzyskania z obwodu z ryc. 1. Tak byłoby, gdyby na przykład założyć, że ucho ma zerową czułość między, powiedzmy, 20 kHz a częstotliwością Nyquista i wybiera funkcję ważenia, aby odzwierciedlić ten fakt.

Ich ryc. 1. pokazuje kształtownik szumów z ogólną strukturą filtra IIR zarówno z biegunami, jak i zerami, tak różny od struktury FIR, którą masz w tej chwili, ale to, co mówią, dotyczy również tego, ponieważ odpowiedź impulsowa filtra FIR może być wykonane dowolnie blisko odpowiedzi impulsowej dowolnego stabilnego filtra IIR.

Skrypt Octave do projektowania filtrów

Oto skrypt Octave do obliczania współczynników inną metodą, która moim zdaniem jest odpowiednikiem Stojanovici i in. Metoda 2014 sparametryzowana jako z właściwym wyborem mojego parametru.$\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

Zaczyna się od okna Dolpha-Chebysheva jako współczynników, zwołuje je, aby podwoić zera funkcji przenoszenia, dodaje do środkowego stuknięcia liczbę, która „podnosi” odpowiedź częstotliwościową (biorąc pod uwagę, że środkowe stuknięcie jest zerowe), więc że jest wszędzie dodatni, znajduje zera, usuwa zera poza okręgiem jednostki, konwertuje zera na współczynniki (wiodący współczynnik od polyzawsze wynosi 1) i odwraca znak co drugi współczynnik, aby filtr górnoprzepustowy . Wyniki (starszej, ale prawie równoważnej wersji) skryptu wyglądają obiecująco:

Rysunek 9. Odpowiedź częstotliwościowa wielkości filtra z (starszej, ale prawie równoważnej wersji) powyższego skryptu.

Rysunek 10. Wykres bieguna zerowego filtra z (starszej, ale prawie równoważnej wersji) powyższego skryptu.

Współczynniki od (starszy ale prawie równoważne wersja) Powyższy skrypt w notacji: {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}. Liczby są duże, co może prowadzić do problemów numerycznych.

Implementacja kształtowania hałasu w oktawach

Wreszcie zrobiłem własną implementację kształtowania hałasu w Octave i nie mam problemów, jak ty. W oparciu o naszą dyskusję w komentarzach, myślę, że ograniczenie w twojej implementacji polegało na tym, że widmo hałasu zostało ocenione za pomocą prostokątnego okna znanego również jako „brak okienkowania”, które rozdzieliło widmo wysokich częstotliwości na niskie częstotliwości.

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');

Rysunek 11. Analiza spektralna szumu kwantyzacji z powyższej implementacji kształtowania szumu dla oktawy dla stałego zerowego sygnału wejściowego. Oś pozioma: znormalizowana częstotliwość. Czarny: brak kształtowania szumów ( c = [1];), czerwony: oryginalny filtr, niebieski: filtr z sekcji „Skrypt oktawy do projektowania filtrów”.

Rysunek 12. Wyjście w dziedzinie czasu z powyższej implementacji kształtowania szumu dla oktawy dla stałego zerowego sygnału wejściowego. Oś pozioma: numer próbki, oś pionowa: wartość próbki. Czerwony: oryginalny filtr, niebieski: filtr z sekcji „Skrypt oktawy do projektowania filtrów”.

Bardziej ekstremalny filtr kształtujący szumy (niebieski) daje bardzo duże skwantowane wyjściowe wartości próbek nawet przy zerowym wejściu.