О погрешности конвольвера

[_Сам_]

Заинтересовала меня недавняя дискуссия в теме «О фазолинейности. Продолжение или революция?!» о погрешности, вносимой конвольвером. После слов ViktKors

Увы, если я слышу грязь, я интерпертирую это как грязь.
Причем, что характрено, если в наушниках я слушаю после конвольвера _искривленную ФЧХ_ , и именно на искривление можно попытаться списать эффект, то в АС ФЧХ ведь выравнивается, а эффект точно так-же плохой.

я решил проверить, какая же «грязь» может возникать после применения фазового FIR-фильтра.

Пусть у нас есть дискретные отсчёты сигнала с частотой выборки 48 кГц. Построим all-pass фильтр с коротким окном на 16 отсчётов. Постараемся сделать так, чтобы коэффициент передачи по амплитуде был равен 1, а фаза менялась по закону косинуса от частоты. При 0 Гц фаза равна 0°, при 12 кГц – минус 90°, а при 24 кГц достигает минус 180°. Для вычисления 16 искомых коэффициентов конвольвера нам нужно решить систему из 16 линейных. Я зафиксировал 7 точек на шкале частот 3,6,9,12,15,18,21 кГц. Для них по условиям задана фаза и амплитуда, что даёт 14 уравнений. Ещё два уравнения даёт условие для коэффиициента передачи +1 на частоте 0 Гц и минус 1 на частоте 24 кГц. Эти частоты стоят особняком, так как на них фаза может быть либо 0° либо ±180°.

Решаем систему и получаем 16 коэффициентов конвольвера:

  0.0000000000000002
  0.9722096263235436
 -0.0000000000000001
 -0.1668931209924746
  0.0000000000000001
 -0.0025951830757338
 -0.0000000000000001
 -0.0016313244799198
  0.0000000000000001
  0.0015642123336106
 -0.0000000000000001
  0.0070858939842292
  0.0000000000000000
  0.0288213444185797
 -0.0000000000000002
  0.1614385514881652

Проверяем, строим график ФЧХ и радуемся:

Целевая фаза показана зелёным цветом, а точки аппроксимации синие кружочки.
Но есть один нюанс. Исходный сигнал содержит, конечно, больше частот, чем наши 9 штук прибитых гвоздями с шагом 3 кГц. Эти частоты тоже буду преобразовываться при свёртке с нашим окном из 16 коэффициентов. Их ФЧХ показана красной кривой. Между узлами интерполяции меняется даже знак производной, т.е. часть частот с положительным ГВЗ будут отставать, а другие с отрицательным ГВЗ – наоборот опережать соседей.

На графике АЧХ картина столь же ужасна.

Хотя точки с шагом 3 кГц точно лежат на целевой зелёной полке 0 дБ, но зато между ними мы огребаем погрешность от -4 до +2 дБ.

Первая мысль была: надо взять больше точек и всё улучшится, но расчёты говорят об обратном.
Увеличил в 16 раз длину конвольвера до 256 точек.


Смотрим и видим, что хотя в самих узлах аппроксимации АФЧХ всё волшебно, но между ними творится кошмар. Одни частоты усиливаются, другие ослабляются, а фаза у всех разъехалась в разные стороны от целевой кривой. Причём амплитуда погрешности и не думает уменьшаться с ростом окна.

Есть у меня подозрения, что именно "лохматость" АФЧХ конвольвера воспринимается на слух как "грязь".

[SilentS]

А как были получены значения АЧХ и ФЧХ в точках, отличных от дискретных 16-и?

[E.Sokol]

Построим all-pass фильтр с коротким окном на 16 отсчётов.

Начнём с того, что это невозможно по определению: all-pass фильтр имеет бесконечную импульсную характеристику, что линейный, что минимально-фазовый каузальный или нет.

---------- Сообщение добавлено 14:11 ---------- Предыдущее сообщение было 14:08 ----------

Её можно обрезать (и наложить гладкое окно при желании) до некоторого кол-ва отсчётов для обеспечения некоторой точности, но 16 для этого явно недостаточно.

[_Сам_]

А как были получены значения АЧХ и ФЧХ в точках, отличных от дискретных 16-и?

дискретная свёртка это просто сумма с весом 16 отстоящих друг от друга по времени на 1/48000 cекунды значений входной функции Fin.
Берём Fin(t) = sin(w * t).
Fout(t) есть сумма 16 слагаемых вида Ci * sin(w*(t-i/48000)) для i=0,15
После собираем амплитуды при sin и cos и вычисляем амплитуду и фазу Fout как функции циклической частоты w.

Начнём с того, что это невозможно по определению: all-pass фильтр имеет бесконечную

пусть невозможно сделать точно. Я моделирую ситуацию с конвольвером, который меняет фазу сигнала по возможности не трогая амплитуду. Ровно это делают при построении фазолинейной АС.

[E.Sokol]

Ровно это делают при построении фазолинейной АС.

Да. Но там берут не 16 отсчётов, а 16К отсчётов, и пульсации АЧХ маскируются шумом квантования выходного устройства, так как ниже его уровня.

---------- Сообщение добавлено 14:42 ---------- Предыдущее сообщение было 14:39 ----------

Аналоговые фильтры типа Чебышева проектируются аналогичным образом - сначала определяется допустимый уровень пульсаций, а уже потом считается сам фильтр.

[_Сам_]

а потому и привел ДВА расчёта. Пульсации АЧХ и при 16 и 256 отсчётах равны 6 дБ. Интересно каким шумом квантования можно это замаскировать?

[E.Sokol]

Интересно каким шумом квантования можно это замаскировать?

В данном случае - однобитным.

---------- Сообщение добавлено 15:21 ---------- Предыдущее сообщение было 14:46 ----------

Я зафиксировал 7 точек на шкале частот 3,6,9,12,15,18,21 кГц.

Это также означает, что амплитуды на всех остальных частотах у вас равны нулю - не удивительно, что allpass фильтр не получился.

---------- Сообщение добавлено 15:23 ---------- Предыдущее сообщение было 15:21 ----------

А пульсации получились вследствие свёртки спектра вашего фильтра со спектром прямоугольной функции.

[_Сам_]

Это также означает, что амплитуды на всех остальных частотах у вас равны нулю - не удивительно, что allpass фильтр не получился.

Дело в том, что у конвольвера с окном в 16 сэмплов и есть всего 16 степеней свободы. Потому всё на что он способен – это контролировать амплитуду и фазу 8 разных частот. Если период входного сигнала тоже составляет 16 выборок, то этого хватает чтобы держать в узде ВСЕ его 8 гармоник. Однако если входной сигнал имеет период 256 отсчётов, то там появятся частоты между 0 и 3 кГц, у которых будет всё плохо с фазой и амплитудой. Суть в том, что даже окна 16к получается мало. Возьмём, к примеру, музыкальный фрагмент длиной 3 минуты и закольцуем его. Получим сигнал с периодом 3*60*48000 = 8.64 млн. сэмплов. Он однозначно преобразуется в 4.32 млн. гармоник от 0 до 24 кГц. Так вот конвольвер длиной 16к, может точно рулить всего 8к гармониками, а остальные 4.312 млн. попадут между узлами и будут линейно искажены.

[E.Sokol]

Дело в том, что у конвольвера с окном в 16 сэмплов и есть всего 16 степеней свободы. Потому всё на что он способен – это контролировать амплитуду и фазу 8 разных частот.

Он не контролирует ни фазу, ни амплитуду. Он производит свёртку с дискретным сигналом у которого все отсчёты, кроме 16 равны нулю.

Однако если входной сигнал имеет период 256 отсчётов, то там появятся частоты между 0 и 3 кГц, у которых будет всё плохо с фазой и амплитудой.

В математике нет понятия "плохо". Значения амплитуды и фазы в этих частотах строго определены. Настолько строго, что их можно посчитать аналитически и они будут совпадать с замеренными.

Суть в том, что даже окна 16к получается мало. Возьмём, к примеру, музыкальный фрагмент длиной 3 минуты и закольцуем его.

Не надо ничего кольцевать. Свёртка в конвольвере - линейная.

[_Сам_]

Он производит свёртку с дискретным сигналом у которого все отсчёты, кроме 16 равны нулю.

а мне удобнее представлять, что это у окна конвольвера кроме 16 коэффициентов все остальные равны нулю, а сигнал может быть любым, даже апериодическим. Окно скользит по оси времени.

В математике нет понятия "плохо".

ну вы поняли о чём я. Плохо значит велико отклонение от целевой функции. Я посчитал сейчас для 1024 выборок. Вся та же ботва. Нет сходимости к allpass. Отклонения не уменьшаются. 6 дБ на АЧХ остаются. Причём фигово, что это близко лежащие частоты. Не ожидал такого. Может я выбрал на редкость неудачный закон изменения фазы.

Не надо ничего кольцевать. Свёртка в конвольвере - линейная.

Закольцовывать надо. Иначе ДПФ не применить.

[E.Sokol]

Первая мысль была: надо взять больше точек и всё улучшится, но расчёты говорят об обратном.

Всё действительно улучшится, если расчёты делать корректно.

Итак, берём 256 отсчётов и рисуем желаемые АЧХ и ФЧХ:


Делаем обратное преобразование Фурье и получаем импульс:


Он (естественно) закольцован и для линейной свёртки не пригоден. Будем раскольцовывать. Сначала сместим (циклически!) к центру:


Затем выравняем левый край с нулём (путём добавления постоянной составляющей) и наложим оконную функцию - это и будет наш FIR-фильтр:


Посмотрим спектр после наложения окна:


Вроде норм, но радоваться рано - поскольку мы ещё в цикле. А вот чтобы увидеть промежуточные частоты - нужно добить FIR-фильтр нулями по краям. В теории надо до бесконечности, на практике достаточно в несколько раз больше исходного размера - добьём до 16К:


Смотрим спектр:


АЧХ не совсем ровная, просто этого не видно. Чтобы увидеть больше, посмотрим разницу между замеренной АЧХ и принудительно выставленной в 0 дБ:


Как видим, никаких пульсаций особо не наблюдается, а максимальное отклонение в -60 дБ. Отклонение ФЧХ от желаемой можно померить аналогичным образом.

---------- Сообщение добавлено 16:53 ---------- Предыдущее сообщение было 16:51 ----------

Эти же манипуляции можно сделать аналитически - будет куча формул, никому не интересных.

[_Сам_]

а можете прикрепить вектор с 256ю коэффициентами окна?

[E.Sokol]

пожалуйста

0
8.0670003004334E-12
1.16561715423389E-12
3.0140270947661E-11
3.14338906155009E-12
7.10512203946516E-11
6.81765708840962E-12
1.45440605243607E-10
1.32174540842347E-11
2.71158902688186E-10
2.35901288241328E-11
4.70689505798514E-10
3.9476030928958E-11
7.72716018335736E-10
6.27851164606949E-11
1.21388016265279E-09
9.58750374718089E-11
1.84078251236361E-09
1.41631644277517E-10
2.71228367968379E-09
2.03551136169692E-10
3.90217220414957E-09
2.85824540768718E-10
5.50227507344301E-09
3.93423719238431E-10
7.62609894833289E-09
5.32189548629675E-10
1.04131077101343E-08
7.0892005044839E-10
1.4033757727544E-08
9.31461062453461E-10
1.86954429644119E-08
1.20879522226136E-09
2.46495255828624E-08
1.5511310803896E-09
3.21996700704753E-08
1.96999255608382E-09
4.17117492085162E-08
2.47830571511484E-09
5.3625643749074E-08
3.09048059847271E-09
6.84693450157091E-08
3.82249948591481E-09
8.68758565297695E-08
4.6919898321683E-09
1.09603522699116E-07
5.7183057308152E-09
1.37560577623405E-07
6.92260285780087E-09
1.71834901648794E-07
8.32790866348122E-09
2.13730249074153E-07
9.95920082427073E-09
2.64810577769733E-07
1.18434732558491E-08
3.26954532847035E-07
1.40098250539957E-08
4.02422801147237E-07
1.64895363283945E-08
4.93941823484932E-07
1.93161708998147E-08
6.04808453227426E-07
2.25256769043986E-08
7.39021653310548E-07
2.61565386723305E-08
9.01449306858573E-07
3.02499110285863E-08
1.09804097194251E-06
3.48498373283452E-08
1.33610132365234E-06
4.00034826769882E-08
1.62464434851028E-06
4.57614296738667E-08
1.97485601746698E-06
5.21780417413358E-08
2.40070416749893E-06
5.93119032839349E-08
2.9197501180729E-06
6.72263145432582E-08
3.55424007306895E-06
7.59899059148709E-08
4.33258894495084E-06
8.56772778144281E-08
5.29142218379432E-06
9.63696589106308E-08
6.47842201820629E-06
1.0815551324574E-07
7.95635124656975E-06
1.21130790071704E-07
9.80883008055575E-06
1.35398595433966E-07
1.21487705006626E-05
1.51067419415419E-07
1.513092124999E-05
1.68246546963892E-07
1.89709137056821E-05
1.87036104094214E-07
2.39748448522975E-05
2.07506438253766E-07
3.05864159272513E-05
2.29656679454236E-07
3.94642393912315E-05
2.53331195013885E-07
5.16128357300011E-05
2.78049244091524E-07
6.86130843776139E-05
3.0264779522555E-07
9.30456108839636E-05
3.24503084196714E-07
0.000129309198548835
3.37744662651172E-07
0.000185301267422728
3.2887415468956E-07
0.00027612960785762
2.65056186869253E-07
0.000433081598408092
5.91635443399624E-08
0.000727898283259952
-5.50797346968178E-07
0.0013488278157209
-2.5045734189668E-06
0.00288952770034781
-1.00906622706341E-05
0.00775675896860692
-4.94218561995732E-05
0.0292615244704905
-0.000317268915708142
0.162289947470143
-0.00195312490616573
0.971484398695247
0.00274493536187026
-0.167564533265104
-0.000399308467979331
-0.00318208863821
-2.07860458959659E-06
-0.00294543965677379
-4.98714582173744E-06
-0.00134370903806331
-2.37744354560726E-06
-0.000725925339764466
-1.43688371515894E-06
-0.000432093484131264
-9.90140612867603E-07
-0.000275586086244994
-7.42950740818166E-07
-0.000184980392461629
-5.89205795874427E-07
-0.000129109225791026
-4.84668693837662E-07
-9.29155607902072E-05
-4.08599936227542E-07
-6.85255600183695E-05
-3.50317365623576E-07
-5.15522571883178E-05
-3.0387785735286E-07
-3.94213245921445E-05
-2.6575519453546E-07
-3.05554152547626E-05
-2.33740252234606E-07
-2.39520772280078E-05
-2.06383104572525E-07
-1.89539547631379E-05
-1.82693080662414E-07
-1.51181346612825E-05
-1.61969026260863E-07
-1.21390280094369E-05
-1.43699012405865E-07
-9.80133900740419E-06
-1.27498686460279E-07
-7.95054536884002E-06
-1.13072114923129E-07
-6.47389090876341E-06
-1.00186062692223E-07
-5.28786440678706E-06
-8.86527936482283E-08
-4.32978052436105E-06
-7.83181704081353E-08
-3.55201282734975E-06
-6.90533206861136E-08
-2.91797656841518E-06
-6.07486800146663E-08
-2.39928686361244E-06
-5.3309683741063E-08
-1.97371990426338E-06
-4.66535819046806E-08
-1.62373119957915E-06
-4.07071141928902E-08
-1.3353657072656E-06
-3.54047154320707E-08
-1.09744722887432E-06
-3.06872080154505E-08
-9.00969311507776E-07
-2.65007858655736E-08
-7.3863311380269E-07
-2.27962425494074E-08
-6.04493629443091E-07
-1.95283818705103E-08
-4.93686553672776E-07
-1.66555549910246E-08
-4.02215732199642E-07
-1.41393112368422E-08
-3.26786535955878E-07
-1.19441098799922E-08
-2.64674299930786E-07
-1.00370969567737E-08
-2.13619745831712E-07
-8.3879290006777E-09
-1.71745359388858E-07
-6.96861207076174E-09
-1.37488093915225E-07
-5.75337129945409E-09
-1.09544922184394E-07
-4.71854470087585E-09
-8.68285566616315E-08
-3.84246653456467E-09
-6.84312428326453E-08
-3.10537722047814E-09
-5.35950216284483E-08
-2.48932167820302E-09
-4.16872055249727E-08
-1.97806086352462E-09
-3.21800609030804E-08
-1.55697709605027E-09
-2.46339152776068E-08
-1.21298034511111E-09
-1.86830681467627E-08
-9.34417081660067E-10
-1.40239944923133E-08
-7.10976596532835E-10
-1.04054469292309E-08
-5.33595639410988E-10
-7.620125540744E-09
-3.94366575117298E-10
-5.49765085121428E-09
-2.86442283929339E-10
-3.89862239916194E-09
-2.03945320556232E-10
-2.70958522569078E-09
-1.41875214246755E-10
-1.83875484179891E-09
-9.6019901578302E-11
-1.21237727922313E-09
-6.2867237165518E-11
-7.71620072876401E-10
-3.95198814036306E-11
-4.6990540862389E-10
-2.36117934506575E-11
-2.7060965972596E-10
-1.32271100137087E-11
-1.4506346303749E-10
-6.82136135178195E-12
-7.07941686394274E-11
-3.14452609546198E-12
-2.99591088592116E-11
-1.16582662284044E-12
-7.92254376208708E-12

[свернуть]

Только ФЧХ здесь слегка отличается от вашей, потому что мне было лень точно подбирать формулу.

[_Сам_]

Только ФЧХ здесь слегка отличается от вашей, потому что мне было лень точно подбирать формулу.

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

у меня получается, что каждая вторая точка сдвинута на 180 градусов относительно нормы.
либо мы как-то по-разному сворачиваем сигналы.

---------- Сообщение добавлено 16:03 ---------- Предыдущее сообщение было 15:33 ----------

кажется понял. Это от постоянного сдвига по оси времени

[E.Sokol]

Это от постоянного сдвига по оси времени

Да, нужно циклически повернуть импульс на пол-экрана перед FFT.

[_Сам_]

Я так понял, у вас свёртка делается центрально-симметрично со сдвигом на полсэмпла. Получается 128 выборок до текущего момента и столько же после него.
В этом случае фаза выглядит так:

Фаза крутится нв 270 градусов вместо 180, но не суть важно.
А у меня была свёртка строго назад по времени без подглядывания в будущее сигнала.
Симметричный подход мне больше нравится

[E.Sokol]

свёртка строго назад по времени без подглядывания в будущее сигнала.

В таком случае нужно проектировать IIR, а не FIR. Через передаточную функцию, Z-преобразование и декомпозицию на биквады.

---------- Сообщение добавлено 23:06 ---------- Предыдущее сообщение было 22:06 ----------

Я так понял, у вас свёртка делается центрально-симметрично со сдвигом на полсэмпла. Получается 128 выборок до текущего момента и столько же после него.
В этом случае фаза выглядит так:

Фаза крутится нв 270 градусов вместо 180, но не суть важно.

Нет, фаза выглядит так, как на моих картинках, и на пол-пикселя я тоже ничего не сдвигаю. Сформировать однонаправленный allpass-импульс через обратное FFT вероятно можно, но я такого способа не знаю. Поэтому он посередине.

[_Сам_]

на пол-пикселя я тоже ничего не сдвигаю.

да, я чуть позже сам понял, что перемудрил. Мой дополнительный сдвиг на полсэмпла как раз внёс линейную составляющую в график фазы (постоянная ГВЗ), которая на 24 кГц достигла дополнительных 90°.

У меня получилось, что отклонение фазы вашего allpass от заданной кривой не превосходит 0.5°


Сейчас попробую своим методом посчитать для симметричной свёртки

---------- Сообщение добавлено 10:58 ---------- Предыдущее сообщение было 09:12 ----------

Симметрия при свёртке рулит.
Сделал длину окна нечётнымм числом, чтобы у него была центральная точка.
Даже при N=15 (Всего 15, Карл!) решение задачи из 1-го поста получается весьма точным.

Вычисленный вектор конвольвера

0.0029219093990154
0.0000000000000000
0.0077736116277327
-0.0000000000000000
0.0291556351355309
0.0000000000000000
0.1615393569250373
-0.0000000000000000
0.9721088301119815
0.0000000000000001
-0.1672274024974516
-0.0000000000000001
-0.0032828881559063
-0.0000000000000000
-0.0029896014495861

[свернуть]

Ошибка амплитуды не превосходит 0.02 дБ

А ошибка фазы 0.4°

Самое приятное, что появилась сходимость. При окне 255 сэмплов ошибки ожидаемо уменьшаются.

[E.Sokol]

Даже при N=15 (Всего 15, Карл!) решение задачи из 1-го поста получается весьма точным.

Секрет в синусе при описании ФЧХ. Возьмите какую-нибудь другую кривую, а лучше ломанную - и результат будет уже не таким обнадёживающим.

[Игорь Гапонов]

Начнём с того, что это невозможно по определению: all-pass фильтр имеет бесконечную импульсную характеристику, что линейный, что минимально-фазовый каузальный или нет.

Правильно зачёркнуто Вообще, для наших баранов попытка шото объяснить с позиций расположения корней полиномов - тупой тупик. Это мой жизненный опыт раздумий. Т.е. закоренелое имхо.

Секрет в синусе при описании ФЧХ. Возьмите какую-нибудь другую кривую, а лучше ломанную - и результат будет уже не таким обнадёживающим.

Секрет в сигнале. Для любого окна, отличного от "бесконечного прямоугольного", любого округления и допуска по точности, отличных от нуля, найдётся сигнал-убийца.