Я хочу сделать формирование шума в 16-битном приложении 100 кГц, чтобы сместить весь шум квантования в полосу 25 кГц-50 кГц с минимальным шумом в полосе постоянного тока 25 кГц.

Я настроил mathematica для создания ядра фильтра ошибок из 31 выборки с помощью обучения с подкреплением, которое хорошо работает: после небольшого обучения я могу получить усиление высокочастотного шума на ~ 16 дБ примерно для того же уровня снижения в низкочастотной полосе ( центральная линия - уровень шума дизеринга). Это соответствует теореме о формировании шума Герцона-Крейвена.

Теперь к моей проблеме:

Мне не удается сформировать шум еще больше даже после обширного обучения, хотя теорема Герцона-Крейвена не запрещает этого. Например, должно быть возможным добиться снижения на 40 дБ в нижней полосе и усиления на 40 дБ в верхней полосе.

Так есть ли еще один фундаментальный предел, с которым я сталкиваюсь?

Я попытался изучить теоремы Шеннона о шуме / дискретизации / информации, но, поиграв с этим некоторое время, мне удалось вывести из нее только один предел: теорему Герзона-Крейвена, которая, кажется, является прямым следствием теоремы Шеннона.

Любая помощь приветствуется.

РЕДАКТИРОВАТЬ: больше информации

Прежде всего, ядро ​​фильтра, которое производит вышеупомянутое формирование шума, Обратите внимание, что самая последняя выборка находится на правой стороне. Числовые значения BarChart округлены до 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} (не совсем гистограмма, но выдает аналогичную кривую )

Еще одно замечание о реализации отзывов об ошибках:

Я попробовал две разные реализации обратной связи об ошибках. Сначала я сравнил округленную выходную выборку с требуемым значением и использовал это отклонение как ошибку. Во-вторых, я сравнил округленную выходную выборку с (вход + ошибка обратной связи). Хотя оба метода дают совершенно разные ядра, оба, похоже, выравниваются примерно с одинаковой интенсивностью формирования шума. Данные, размещенные здесь, используют вторую реализацию.

Вот код, который используется для расчета оцифрованных образцов волны. шаг - размер шага для округления. волна - это нецифрованный сигнал (обычно это нули, когда сигнал не подается).

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
  ]

Метод армирования:

«Балл» рассчитывается путем анализа спектра мощности шума. Цель состоит в том, чтобы минимизировать мощность шума в диапазоне DC-25 кГц. Я не наказываю шум в высокочастотном диапазоне, поэтому произвольно высокий уровень шума не уменьшит оценку. Я ввожу шум в веса ядра, чтобы учиться. Возможно, поэтому я нахожусь в (очень широком и глубоком) локальном минимуме, но я считаю это крайне маловероятным.

Сравнение со стандартным дизайном фильтра:

Mathematica позволяет создавать фильтры итеративно. Они могут иметь намного лучший контраст, чем 36 дБ, когда их частотная характеристика построена; до 80-100 дБ. Числовые значения: {0,024, -0,061, -0,048, 0,38, -0,36, -0,880, 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,880, -0,36, 0,38, -0,048, -0,061, 0,024}

Тем не менее, когда они применяются при фактическом формировании шума, они (а) фиксируются с одинаковым контрастом ~ 40 дБ, (б) работают хуже, чем запомненный фильтр, фактически не ослабляя шум.

Основное дизеринг без шумоподавления

Базовое дитрированное квантование без ограничения шума работает следующим образом:

Рисунок 1. Базовая диаграмма системы дитризованного квантования. Шум - это треугольное сглаживание с нулевым средним с максимальным абсолютным значением 1. Округление до ближайшего целого числа. Остаточная ошибка - это разница между выходом и входом, она рассчитывается только для анализа.

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

С независимой аддитивной остаточной ошибкой мы получили бы более простую модель системы:

Рисунок 2. Аппроксимация основного дитризованного квантования. Остаточная ошибка - белый шум.

В приближенной модели выходной сигнал является просто входом плюс независимая остаточная ошибка белого шума.

Дизеринг с формированием шума

Я не могу читать Mathematica очень хорошо, поэтому вместо вашей системы я проанализирую систему из Lipshitz et al. « Минимально слышимое формирование шума » J. Audio Eng. Soc., Vol.39, No.11, November 1991:

Рисунок 3. Липшиц и соавт. 1991 системная схема (адаптировано из их рис. 1). Фильтр (выделенный курсивом в тексте) включает в себя задержку в одну выборку, чтобы его можно было использовать в качестве фильтра обратной связи по ошибкам. Шум треугольный дизеринг.

Если остаточная ошибка не зависит от текущих и прошлых значений сигнала A, мы имеем более простую систему:

Рисунок 4. Примерная модель Липшица и соавт. Система 1991 года. Фильтр такой же, как на рис. 3, и включает в себя задержку в одну выборку. Он больше не используется в качестве фильтра обратной связи. Остаточная ошибка - белый шум.

В этом ответе я буду работать с более легко анализируемой приближенной моделью (рис. 4). В оригинале Lipshitz et al. Система 1991 года, Filter имеет общую форму фильтра с бесконечной импульсной характеристикой (IIR), которая охватывает как фильтры IIR, так и фильтры с конечной импульсной характеристикой (FIR). В дальнейшем мы будем предполагать, что Фильтр является КИХ-фильтром, поскольку, как я полагаю, исходя из моих экспериментов с вашими коэффициентами, именно это вы имеете в своей системе. Передаточная функция фильтра :

$z^{-1}$

$\ldots, -b_3, -b_2, -b_1$$1, b_1, b_2, b_3, \ldots$$b_0 = 1$horzcatв скрипте Octave ниже), и, наконец, список переворачивается (путем 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)

Скрипт отображает амплитудно-частотную характеристику и нулевые местоположения фильтра полного шумоподавления:

Рисунок 5. АЧХ полного шумоподавляющего фильтра.

$\times$$\circ$

Я думаю, что проблема нахождения коэффициентов фильтра может быть переформулирована как проблема разработки фильтра минимальной фазы с ведущим коэффициентом 1. Если существуют частотные характеристики таких фильтров, то эти ограничения переносятся в эквивалентные ограничения. в формировании шума, который использует такие фильтры.

Переход от всеполюсной конструкции к минимально-фазной РПИ

Процедура проектирования различных, но во многих отношениях эквивалентных фильтров описана в Stojanović et al. , "Проектирование всеполюсных рекурсивных цифровых фильтров на основе сверхсферических полиномов", Radioengineering, том 23, № 3, сентябрь 2014 года. Они рассчитывают знаменательные коэффициенты передаточной функции всеполюсного фильтра низких частот IIR. Они всегда имеют ведущий коэффициент знаменателя, равный 1, и имеют все полюса внутри единичного круга, что является требованием стабильных БИХ-фильтров. Если эти коэффициенты используются в качестве коэффициентов фильтра формирования шума КИХ с минимальной фазой, они будут давать инвертированную частотную характеристику верхних частот по сравнению с фильтром БИХ нижних частот (знаменательные коэффициенты передаточной функции становятся коэффициентами числителя). В вашей записи есть один набор коэффициентов из этой статьи {-0.0076120, 0.0960380, -0.5454670, 1.8298040, -3.9884220, 5.8308660, -5.6495140, 3.3816780}, который можно протестировать для приложения формирования шума, хотя это не совсем соответствует спецификации:

Рисунок 7. АЧХ фильтра FIR с использованием коэффициентов от Стояновича и соавт. 2014.

Рисунок 8. График FIR-фильтра с нулевым полюсом с использованием коэффициентов Стояновича и соавт. 2014.

Всеполюсная передаточная функция:

$a$$b$

Чтобы спроектировать всеполюсный фильтр и преобразовать его в КИХ-фильтр с минимальной фазой, вы не сможете использовать методы расчета БИХ-фильтра, которые начинаются с аналогового фильтра-прототипа и отображают полюсы и нули в цифровую область с помощью билинейного преобразования. , Это включает в себя cheby1, cheby2и ellipв Октаве и SciPy Python. Эти методы будут отдавать нули от начала координат z-плоскости, поэтому фильтр не будет иметь требуемого всеполюсного типа.

Ответ на теоретический вопрос

Если вам все равно, сколько шума будет на частотах выше четверти частоты дискретизации, тогда Lipshitz et al. 1991 отвечает на ваш вопрос напрямую:

Для таких весовых функций, которые стремятся к нулю в части полосы, теоретическое ограничение для взвешенного снижения мощности шума, получаемое по схеме на рис. 1, отсутствует. Это было бы в том случае, если, например, предполагается, что Ухо имеет нулевую чувствительность, скажем, между 20 кГц и частотой Найквиста, и выбирает весовую функцию для отражения этого факта.

На их Рис. 1 показан формирователь шума с общей структурой фильтра БИХ с обоими полюсами и нулями, настолько отличающийся от структуры КИХ, которая у вас есть в данный момент, но то, что они говорят, относится и к этому, потому что импульсная характеристика фильтра КИХ может быть сделано произвольно близко к импульсной характеристике любого данного стабильного БИХ-фильтра.

Октавный скрипт для дизайна фильтров

$\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

Он начинается с окна Дольфа-Чебышева в качестве коэффициентов, сворачивает его с собой, чтобы удвоить нули передаточной функции, добавляет к среднему отводу число, которое «поднимает» частотную характеристику (рассматривая среднее отвод как нулевое время), так что он везде положительный, находит нули, удаляет нули, находящиеся за пределами единичного круга, преобразует нули обратно в коэффициенты (начальный коэффициент polyвсегда равен 1) и отображает знак каждого второго коэффициента, чтобы сделать фильтр высокочастотным , Результаты (более старой, но почти эквивалентной версии) сценария выглядят многообещающе:

Рисунок 9. Амплитудно-частотная характеристика фильтра из (более старой, но почти эквивалентной версии) приведенного выше сценария.

Рис. 10. График фильтра с нулевого полюса из (более старой, но почти эквивалентной версии) приведенного выше сценария.

Коэффициенты из (старая , но почти эквивалентна версия) выше сценария в обозначениях: {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}. Числа велики, что может привести к численным проблемам.

Октавная реализация шумоподавления

Наконец, я сделал свою собственную реализацию шумоподавления в Octave, и у меня не возникло таких проблем, как у вас. Основываясь на нашем обсуждении в комментариях, я думаю, что ограничение в вашей реализации заключалось в том, что спектр шума оценивался с использованием прямоугольного окна, называемого «без окон», которое проливало высокочастотный спектр на низкие частоты.

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

Рис. 11. Спектральный анализ квантования шума из описанной выше октавной реализации формирования шума для входного сигнала с постоянным нулем. Горизонтальная ось: нормализованная частота. Черный: нет формирования шума ( c = [1];), красный: ваш оригинальный фильтр, синий: фильтр из раздела «Скрипт октавы для дизайна фильтра».

Рисунок 12. Выход во временной области из описанной выше реализации октавного формирования шума для входного сигнала с постоянным нулем. Горизонтальная ось: номер образца, вертикальная ось: значение образца. Красный: ваш оригинальный фильтр, синий: фильтр из раздела «Октавный скрипт для дизайна фильтра».

Фильтр с более экстремальным шумоподавлением (синий) приводит к очень большим квантованным значениям выходной выборки даже при нулевом входном сигнале.