Dithering (часть 3)

Это будет третья и завершающая часть из серии постов про dithering, да и тема, наверняка, уже порядком надоела ;)

Шум для dithering'а, о котором я говорил в предыдущих постах был либо белый, либо топорно приглушенный в нижних частотах. Технология noise shaping, о которой можно прочитать в статье S.P. Lipshitz, J. Vanderkooy and R.A. Wannamaker. 1991. Minimally audible noise shaping. J. Audio Eng. Soc. 39, 836-852. — это работа скальпеля, а не топора. Дело в том, что человек воспринимает разные частоты неодинаково. Существуют экспериментальные кривые зависимости мощности тона от частоты, которые построены так, что человек воспринимает эти тоны с одинаковой "громкостью". Вот пример такой кривой (спасибо, Википедия ;-)

Как видно, на частотах около 4 и 12 кГц у человека есть локальные "пики" восприимчивости, когда тон должен быть тише, чем тоны с соседними частотами, чтобы восприниматься субъективно с той же громкостью.

Схема dithering'а выглядит так (взял из упомянутой статьи):

К исходному сигналу добавляется белый шум (например, треугольно распределенный) и сигнал квантуется. Затем из полученного зашумленно-квантованного сигнала вычитается исходный (находится ошибка округления), она фильтруется, так, чтобы шум имел желаемый спектр (в некотором смысле он должен быть обратной величиной к кривой Perceived Human Hearing выше, чтобы человек воспринимал все частотные составляющие одинаково). Если я правильно понял фишку (могу и ошибаться), то при подборе фильтра учитывается ещё то, что минимизируется "взвешенная" мощность ошибки+шума $W_{w} = \int_{-\omega_{ref}/2}^{\omega_{ref}/2} W(\omega) A(\omega) d\omega$, где $W(\omega)$ — спектральная плотность мощности, $A(\omega)$ — её вес с учетом человеческого восприятия.

В частности, в программе audacity применяется фильтр ошибки, ядро которого имеет пять ненулевых коэффициентов (скопипащенных из той же статьи, копипаст — великая сила) $a_{n} = \left\{2.033, -2.165, 1.959, -1.590, 0.6149, 0, 0, \cdots \right\}$ Ошибка, вычисленная таким образом для $n$-го семпла, складывается потом с входным сигналом на $n+1$ шаге. Вот код, спертый из audacity:

(defun uniform-noise (magnitude)
  (- (random (ash 1 magnitude))
     (ash 1 (1- magnitude))))

(defun triangle-noise (magnitude)
  (+ (uniform-noise magnitude)
     (uniform-noise magnitude)))

(defun apply-dither (array magnitude noise-fn)
  (let ((new-array (make-array (length array))))
    (map-into new-array
              (lambda (x)  (funcall noise-fn magnitude x))
              array)))

(defun quantize (x m)
  (let ((rnd (ash 1 m)))
    (* rnd (round x rnd))))

(defun shaped-noise ()
  (let ((history (make-array 8 :initial-element 0))
        (flt #(2.033 -2.165 1.959 -1.590 0.6149))
        (position 0))
    (lambda (magnitude sample)
      (let* ((noise (triangle-noise magnitude))
             (xe (+ sample
                    (* (aref flt 0) (aref history (logand position 7)))
                    (* (aref flt 1) (aref history (logand (- position 1) 7)))
                    (* (aref flt 2) (aref history (logand (- position 2) 7)))
                    (* (aref flt 3) (aref history (logand (- position 3) 7)))
                    (* (aref flt 4) (aref history (logand (- position 4) 7)))))
             (result (+ xe noise)))
        (prog1 result
          (setq position (logand (1+ position) 7))
          (setf (aref history position)
                (- xe (quantize result magnitude))))))))

Таким образом, сам шум получается отфильтрованным эквивалентным фильтром с ядром $\left\{-1, 2.033, -2.165, 1.959, -1.590, 0.6149, 0, 0, \cdots \right\}$. Вот АЧХ этого фильтра (вертикальная шкала — логарифмичекая):

АЧХ эквивалентного фильтра для шума. По горизонтальной оси — частоты в долях $\omega_{ref}$

Взяв частоту дискретизации $f_{ref} = 44100 Hz$, можно подсчитать, что локальные минимумы находятся около частот $f_{1} \approx 0.1 f_{ref} \approx 4000 Hz$ и $f_{2} \approx 0.3 f_{ref} \approx 13000 Hz$. Это идеально согласуется с картинкой в начале поста — раз человек слышит 4 кГц лучше, то имеет смысл немного "приглушить" эту частоту, чтобы она не выделялась на фоне других. Также видим, что верхние частоты звучат сильнее нижних (как в фильтрованном треугольном шуме в предыдущем посте).

У этого метода есть свой недостаток, однако. Раз какая-то часть энергии шума будет сдвинута в верхние частоты за счет нижних, то и амплитуда шума будет больше, чем при белом шуме (для компенсации "приглушенности" нижнечастотного шума). Её можно оценить как $A = 2\Delta (1 + \sum_{k=1}^{5} \left|a_{k}\right| ) \approx 9 \cdot 2\Delta$. В вики audacity рекомендуется применять этот вид dithering'а не более одного раза. Также не совсем ясно, что делать, если значение сигнал+шум превысит допустимые значения. Audacity просто клиппует результирующий сигнал, но не внесет ли это в данном случае бОльшую ошибку, чем та, от которой мы избавляемся (в немного меньшей мере эта проблема относится ко всем видам шума)?

Даю послушать тот же сигнал, что и всегда, но с примененным noise shaping:
Sine wave (shaped dithering, 5 coeffs) 400Hz (44100 Hz, 7 bps)

На спектрограмме видим:

Весь шум здесь сдвинут в верхние частоты, а потому звучит тихо.

Также в статье есть пример фильтра с ядром, имеющим 3 ненулевых коэффициента: $\left\{ 1.652, -1.049, 0.1382, 0, \cdots, \right\}$. Сам фильтр имеет худшие психоакустические характеристики, но и амплитуда шума у него меньше.

Пример:
Sine wave (shaped dithering, 3 coeffs) 440 Hz (44100 Hz, 7 bps)

Спектр:

В заключение хочу напомнить, что все эти коэффициенты имеют смысл только при частоте дискретизации $f_{ref} = 44100 Hz$! При других частотах теряется весь смысл такого фильтра (локальные минимумы на АЧХ эквивалентного фильтра для шума сдвигаются), и коэффициенты надо считать заново. Audacity, однако, не учитывает этого. Поэтому, при сохранении чего-либо в audacity с частотой дискретизации, отличной от 44100 Hz, лучше использовать triangular noise. Сказанное относится к версии 2.2.2