Как правильно реализовать КИХ-эквалайзер? Опции

[elk]

Такой вот вопрос. Я реализую алгоритм свертки и хочу прикрутить к ней дополнительно эквалайзер, причем чтобы форму кривой АЧХ эквалайзера можно было задавать в принципе любой (т.е. не обязательно привязываться к комбинации стандартных фильтров типа Баттерворта и т.д.). Свертка, естественно, реализуется через БПФ (допустим, overlap-add). Входной сигнал поступает на обработку кусками, для простоты ограничимся кусками, равными целой степени двойки (напр., 512, 1024).
Изначальная идея состоит в том, что раз фильтр задается в частотной области, то нужно просто домножить результат преобразования Фурье (ПФ) ИХ на характеристику эквалайзера (т.е. каждый комплексный коэффициент ПФ ИХ на значение коэффициента усиления эквалайзера на данной частоте). Это будет эквивалентно тому, что комплексные коэффициенты БПФ импульсной характеристики эквалайзера являются полностью действительными. Это, по-видимому, не совсем правильно, т.к. фильтры как правило моделируют минимально-фазовыми цепями, у которых АЧХ и ФЧХ связаны преобразованием Гильберта. Допустим, мы реализуем преобразование Гильберта, по АЧХ ищем ФЧХ, по ним определяем правильные действительные и мнимые части ПФ эквалайзера.
Вопрос у меня вот в чем. Для получения корректного результата линейной свертки нужно, чтобы оба сигнала во временной области были дополнены нулями до размерности ПФ, так чтобы суммарная длительность ненулевых фрагментов не превышала длину ПФ, иначе получим циклическую свертку. Но ИХ эквалайзера будет ненулевой по всей длине... получается, для корректного результата надо ПФ эквалайзера преобразовать в ИХ эквалайзера с помощью обратного ПФ, обнулить у нее половину отсчетов и вернуть назад в ПФ (причем АЧХ при этом чуть изменится...), и только после этого множить на ПФ ИХ? Но тут возникает та же проблема - мы теперь эту "эквализованную" ИХ будем еще сворачивать с сигналом, т.е. у нее тоже половина отсчетов должна быть нулевой... Опять ОПФ-обнуление коэффициентов-ПФ? Какая-то дикая схема получается... Ощущение такое, что где-то в рассуждениях допущена принципиальная ошибка, только не могу понять где именно...
В результате реализации алгоритма без учета эффекта циклической свертки (т.е. если не играться с ОПФ-обнулением-ПФ, а просто множить ПФ ИХ на ПФ эквалайзера), на краях кусков выходного сигнала (тех которые длиной типа 512, 1024) иногда наблюдаются разрывы (собственно чтобы как-то объяснить их появление и была разведена вся эта теория )

Уточню одну вещь. Мне интересна реализация именно данного принципа (БПФ-свертка + КИХ-эквалайзер), потому советы использовать что-то другое не принимаются

[bahurin]

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

Сообщение отредактировал bahurin - Jan 4 2011, 14:06

[Alexey Lukin]

bahurin, если эквалайзер делать на основе БПФ без использования в явном виде импульсной характеристики, то эквалайзер не будет линейным и инвариантным к сдвигу.

elk, вы сами себя запутали. Вот у вас импульс длины M и блок сигнала длины N. Выбираете размер FFT K так, что K >= M + N - 1. Вычисляете FFT размера K от импульса длины M и сигнала длины N. Перемножаете спектры, берёте обратное FFT, получаете блок результирующего сигнала длиной M + N - 1.

[elk]

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

Хм... Вообще-то именно с этого все и начиналось:

Изначальная идея состоит в том, что раз фильтр задается в частотной области, то нужно просто домножить результат преобразования Фурье (ПФ) ИХ на характеристику эквалайзера (т.е. каждый комплексный коэффициент ПФ ИХ на значение коэффициента усиления эквалайзера на данной частоте).

Я, правда, как бы ИХ эквализую, а не сам сигнал, чтобы делать это всего один раз при изменении АЧХ эквалайзера, а не с каждым новым куском сигнала...

НО надо делать обработку с перекрытием иначе на стыках будут щелчки.

А можно поподробнее?

bahurin, если эквалайзер делать на основе БПФ без использования в явном виде импульсной характеристики, то эквалайзер не будет линейным и инвариантным к сдвигу.

Т.е. я правильно считаю, что нельзя просто множить на АЧХ эквалайзера, а требуются какие-то дополнительные фокусы? Так как, собственно, при этом надо "использовать в явном виде импульсную характеристику"?

Вот у вас импульс длины M и блок сигнала длины N. Выбираете размер FFT K так, что K >= M + N - 1. Вычисляете FFT размера K от импульса длины M и сигнала длины N. Перемножаете спектры, берёте обратное FFT, получаете блок результирующего сигнала длиной M + N - 1.

Без эквалайзера все именно так и работает. Вопрос: как в эту красивую и стройную схему впихнуть эквалайзер, чтоб он ничего не попортил?
Неужели единственный выход - каждый блок сигнала сворачивать с ИХ эквалайзера?
Т.е. после того, как мы получили "блок результирующего сигнала длиной M + N - 1", теперь дополнить его еще кучей нулей (для FFT размером S >= M + N - 1 + L - 1, где L - длина ИХ эквалайзера), вычислить FFT, умножить на FFT эквалайзера и наконец после обратного FFT получить окончательный результат?

(а еще у меня ИХ тоже на блоки побита... это ж как надо схему переколпаковать да перывыколпаковать )

[Alexey Lukin]

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

Правильно. При умножении двух спектров получается круговая свёртка, а не линейная. Поэтому ваша система будет нелинейна и неинвариантна к сдвигу. Будет ли это заметно на практике — зависит от крутизны АЧХ эквалайзера. Но нормальные эквалайзеры так не делают. Они честно сворачивают со своей импульсной характеристикой.

[elk]

Но нормальные эквалайзеры так не делают. Они честно сворачивают со своей импульсной характеристикой.

Ну так сворачивать со своей импульсной можно и ту импульсную характеристику, с которой собственно свертка... Т.е. берем FFT ИХ, множим на FFT эквалайзера, считаем обратное FFT, получаем эквализованную ИХ, дальше уже к ней применяем схему свертки с сигналом...
Потому возвращаюсь к изначальному вопросу: как получить правильное FFT эквалайзера? Правильна ли следующая схема:
АЧХ длиной К - (преобразование Гильберта) - АЧХ+ФЧХ - комплексный спектр - (обратное FFT) - ИХ эквалайзера длиной К - (обнуление К/2 отсчетов, или наоборот добавление К нулевых) - (возможно, еще домножить на окно) - правильная ИХ эквалайзера - (FFT) - правильное FFT эквалайзера длиной К или 2*К (смотря как нули дописывали)

Хотелось, конечно, прилепить эквалайзер к уже готовым блокам FFT ИХ, которые используются в свертке, и мне что-то кажется, что это тоже теоретически возможно (просто бОльшим количеством нулей дополнять блоки ИХ?)

[bahurin]

bahurin, если эквалайзер делать на основе БПФ без использования в явном виде импульсной характеристики, то эквалайзер не будет линейным и инвариантным к сдвигу.

Господь с вами! Как это не будет линейным? Будет линейным никуда не денется. ФЧХ у этого эквалайзера будет нелинейна, это правда, но ее можно задать линейной при необходимости и умножать на комплексный к-т передачи эквалайзера уже с линейной ФЧХ. Для правильной работы такого эквалайзера необходимо делать обработку с перекрытием. Т.е. Берете кусок длиной N берете его fft умножаете на АЧХ берете ifft. Если все правильно сделано мнимая часть ifft должна быть очень маленькой почти 0. Поэтому берем только реальную часть ifft. Из этой реальной части выбираем центральные N/2 отсчетов, т.е от N/4 до 3N/4 - это результат. После сдвигаем исходный сигнала на N/2, снова берем N отсчетов и снова fft-умножение-ifft и выборка центральных N/2 отсчетов. Т.о. вы сдвигаетесь каждый раз на N/2 отсчетов, но делаете N точечное fft и выбираете из результата ifft тоже N/2 отсчетов, дырок нет. В первом буфере первые N/4 отсчетов в результате нули, а центральные N/2 отсчетов которые вы выбираете из результата идут встык друг с другом.

Сообщение отредактировал bahurin - Jan 5 2011, 14:08

[Alexey Lukin]

bahurin, строго говоря, такая схема не обеспечивает инвариантности к сдвигу, а значит — система не будет эквивалентна свёртке (в некоторых случаях — будет порождать гармоники).
Когда FFT размера N от куска данных размера N умножается на АЧХ фильтра, то возникает круговая свёртка, и результирующий сигнал из N точек портится целиком. Так что, извлекая из него "серединку", вы лишь несколько уменьшите артефакты.

[phantom]

В последнее время, что-то всех в электронную музыку понесло, как видно. - постановка задачи мне до боли напоминает музыкальный конволюционный процессор. Ну так в широкоизвестном софтовом Waves IR-1, при изменении настроек эквалайзера происходит пересчет ИХ, что видно по изменении ее изображения. Не согласиться с Alexey Lukin логика не позволяет, но наверное и "elk:Хотелось, конечно, прилепить эквалайзер к уже готовым блокам FFT ИХ, которые используются в свертке, и мне что-то кажется, что это тоже теоретически возможно " - имеет под собой нечто.

[Serg76]

2 Alexey Lukin Вы не верите, что сигнал можно обрабатывать в частотной области

[Alexey Lukin]

Почему же? Очень даже верю! Просто я призываю к аккуратности. Аккуратный эквалайзер — это свёртка. Свёртку можно вычислять в частотной области. Для этого есть 2 метода: overlap-add и overlap-save. То, что описал bahurin — это ни то, ни другое, и вообще не свёртка (т.к. нет ограничений на АЧХ эквалайзера).
Желание elk повторно использовать имеющиеся FFT — понятно! Но пока что нам нигде подробно не рассказали, какие же FFT уже имеются, с каким шагом по времени они идут, что делают и т.д...

[phantom]

Как я понимаю, у elk-a правильно реализована блочная свертка с ИХ, когда входной сигнал и ИХ представлены по-блочно, с разными длинами блоков и вопрос только в том, как преобразовать (отфильтровать) исходную ИХ с использованием эквалайзера (и можно ли это сделать в принципе). Я думаю, схему блочной свертки переменной длины довольно тяжело представить в подробностях наглядно(графически? текстом?), поскольку там наверное должны быть некие комбинации overlap-add или overlap-save методов. Но надо ли это делать вообще? Лучше подождем elk-a с разъяснениями

[Alexey Lukin]

Как я понимаю, у elk-a правильно реализована блочная свертка с ИХ ... и вопрос только в том, как преобразовать (отфильтровать) исходную ИХ с использованием эквалайзера

Если задача действительно такая, то обрабатываем имеющуюся ИХ эквалайзером (в процессе чего она может несколько удлиниться) — и производим свёртку уже имеющимся алгоритмом (но с новым ядром).

[elk]

Берете кусок длиной N берете его fft умножаете на АЧХ берете ifft. Если все правильно сделано мнимая часть ifft должна быть очень маленькой почти 0. Поэтому берем только реальную часть ifft. Из этой реальной части выбираем центральные N/2 отсчетов, т.е от N/4 до 3N/4 - это результат. После сдвигаем исходный сигнала на N/2, снова берем N отсчетов и снова fft-умножение-ifft и выборка центральных N/2 отсчетов. Т.о. вы сдвигаетесь каждый раз на N/2 отсчетов, но делаете N точечное fft и выбираете из результата ifft тоже N/2 отсчетов, дырок нет. В первом буфере первые N/4 отсчетов в результате нули, а центральные N/2 отсчетов которые вы выбираете из результата идут встык друг с другом.

Соглашусь с Alexey Lukin.
Похоже на неправильное описание метода overlap-save.
Проблема в том, что, для избежания циклической свертки, FFT одного из сигналов должно быть получено от блока, содержащего половину нулей.

Как я понимаю, у elk-a правильно реализована блочная свертка с ИХ, когда входной сигнал и ИХ представлены по-блочно, с разными длинами блоков и вопрос только в том, как преобразовать (отфильтровать) исходную ИХ с использованием эквалайзера (и можно ли это сделать в принципе).

Совершенно верно. Насчет разных длин блоков - ну вы, однако, прям мои мысли читаете

Желание elk повторно использовать имеющиеся FFT — понятно! Но пока что нам нигде подробно не рассказали, какие же FFT уже имеются, с каким шагом по времени они идут, что делают и т.д...

Ок, рассказываю подробнее. Имеются FFT от блоков ИХ разного размера (non-uniform partitioned convolution, описанная у Gardner). Самый маленький размер блока ИХ равен длине блока входного сигнала (который может быть 512, 1024 и т.д.). Для каждого размера блоков используется схема overlap-save, т.е. FFT считаются от блоков ИХ, дополненных нулями, а блоки входного сигнала берутся с перекрыванием.
Насчет шага по времени - не совсем понятно, что имеется в виду... Блоки сигнала приходят на обработку поочередно, общая длина сигнала заранее неизвестна. Допустим, что имеется достаточное время для обработки каждого блока согласно описанной схеме.

Если задача действительно такая, то обрабатываем имеющуюся ИХ эквалайзером (в процессе чего она может несколько удлиниться) — и производим свёртку уже имеющимся алгоритмом (но с новым ядром).

Алексей, хотелось бы все-таки услышать Ваш комментарий по поводу правильности способа получения FFT эквалайзера, описанного в моем предыдущем посте...
Или это какая-то страшная тайна ?

[Alexey Lukin]

Я просто не очень понимаю, в чём тут затруднение. Раз имеется ИХ, разбитая на блоки, то (1) составим блоки ИХ вместе, (2) отфильтруем эквалайзером, (3) снова разобьём на блоки — и вперёд, сворачивать со входным сигналом по Гарднеру. Если ядро длинное, то шаг (2) можно делать через FFT, в т.ч. с использованием уже имеющихся FFT блоков исходной ИХ.

Кстати, обратите внимание, что метод Гарднера запатентован.

[elk]

Я просто не очень понимаю, в чём тут затруднение. Раз имеется ИХ, разбитая на блоки, то (1) составим блоки ИХ вместе, (2) отфильтруем эквалайзером, (3) снова разобьём на блоки — и вперёд, сворачивать со входным сигналом по Гарднеру. Если ядро длинное, то шаг (2) можно делать через FFT, в т.ч. с использованием уже имеющихся FFT блоков исходной ИХ.

Кстати, обратите внимание, что метод Гарднера запатентован.

Ну вот этот этап меня интересует:
(2) отфильтруем эквалайзером
Просто умножать на АЧХ эквалайзера нельзя, прокомментировать правильность алгоритма получения ИХ никто не хочет
Еще раз: берем АЧХ, с помощью преобразования Гильберта получаем ФЧХ, из АЧХ и ФЧХ получаем комплексный спектр, из него через IFFT получаем ИХ, обнуляем у нее половину отсчетов, умножаем на окно, считаем FFT - вот ЭТО теперь используем для FFT-свертки с ядром.
Это правильно или нет? Хочу просто ответ на этот вопрос, хотя бы да или нет
Мне еще встречался вариант без Гильберта, когда АЧХ запихивается в действительную часть, а мнимая ставится 0 (т.е. ФЧХ принимается нулевой)... Так правильно? Вроде при этом реальная АЧХ чуть сильнее искажается...
Кстати, если минимальный размер блока очень маленький (например, 64), использовать имеющиеся блоки будет не совсем хорошо, т.к. придется ограничить длину ИХ эквалайзера 64 отсчетами, т.е. сильно исказить его АЧХ... Это рассуждение правильно?

А насчет метода Гарднера, то вроде как Farina вполне убедительно доказывает, что в Lake DSP Patent (Garcia) запатентован не сам метод свертки, а только оптимальный метод разбивки на блоки с помощью алгоритма Витерби...

[petrov]

Просто умножать на АЧХ эквалайзера нельзя

Можно сделать банк фильтров(FFT + полифазный фильтр) и изменять коэффициенты передачи фильтров без ограничений, то бишь "умножать на АЧХ".

[elk]

petrov, а можно что-то типа ссылки на первоисточники?
А то у меня полифазные фильтры как-то узко ассоциировались с ресамплингом...
Я правильно понимаю, что с помощью полифазного фильтра можно именно короткую ИХ эквалайзера более эффективно прицепить, за счет того что для более длинных блоков не надо добавялть лишнюю кучу нулей?
Мне бы более предпочтительно было наоборот подлиннее, типа там 16тыщ отсчетов хотя бы (где-то кто-то писал, но не помню где и кто , что для качественного эквалайзера местами требуется ИХ длиной около 300мс)

Кстати, ни у кого не завалялась статейка S. Lin and S. K. Mitra, “Overlapped block digital filtering”?

[blackfin]

... а можно что-то типа ссылки на первоисточники?

Vaidyanathan P.P.: tutorial.

[Alexey Lukin]

Ну вот этот этап меня интересует:
(2) отфильтруем эквалайзером

Этот этап реализуется в виде свёртки импульса из шага (1) с ИХ эквалайзера. Причём тут преобразование Гильберта — не понял...

А насчет метода Гарднера, то вроде как Farina вполне убедительно доказывает, что в Lake DSP Patent (Garcia) запатентован не сам метод свертки, а только оптимальный метод разбивки на блоки с помощью алгоритма Витерби...

Если мне не изменяет память, в патенте Lake алгоритм Витерби вообще не упоминается, а патентуется именно неравномерное разбиение. А Витерби-разбиение запатентовано в другом патенте, от Microsoft.

Можно сделать банк фильтров(FFT + полифазный фильтр) и изменять коэффициенты передачи фильтров без ограничений, то бишь "умножать на АЧХ".

И каким же образом это поможет решить поставленную задачу? (объединение уже имеющейся секционной свёртки с новым эквалайзером)
Вот и я не знаю.

[petrov]

И каким же образом это поможет решить поставленную задачу? (объединение уже имеющейся секционной свёртки с новым эквалайзером)
Вот и я не знаю.

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

[phantom]

Этот этап реализуется в виде свёртки импульса из шага (1) с ИХ эквалайзера. Причём тут преобразование Гильберта — не понял...

По-моему у elk-a просто затруднения с тем, как получить ИХ по АЧХ эквалайзера, и надо ли при этом учитывать его ФЧХ. Если нет ФЧХ, ее иногда можно получить с помощью преобразования Гильберта. Вот тут http://www.dspguide.com/ch17/1.htm есть вариант, но без учета ФЧХ.

[ivan219]

Сделать ИХ по АЧХ очень просто. Берём все Im = 0 а из Re троим АЧХ после чего делаем обратное FFT получаем временной сигнал сдвигаем его на N/2 отсчётов в право после чего обнуляем 0 отсчёт и накладываем окно. В итоге получается ИХ длинной N-1 отсчёт. Можно также и обрезать ИХ до нужного размера.

В приложении программа и исходники на Delphi может пригодится без учёта ФЧХ.

Сообщение отредактировал ivan219 - Jan 9 2011, 14:06

Прикрепленные файлы

 ______________________________________.rar ( 169.55 килобайт ) Кол-во скачиваний: 66

 

[Alexey Lukin]

Тогда свёртка тут непричём, а задача называется "расчёт фильтра" (filter design). Простейший метод — это действительно метод взвешивания, описанный ivan219.

Кстати, для звукового эквалайзера во многих случаях выгодно использовать биквадратные IIR-фильтры, а не линейнофазовый FIR, получающийся методом взвешивания. Вычислительные затраты обычно получаются значительно меньше.

[phantom]

Тогда свёртка тут непричём, а задача называется "расчёт фильтра" (filter design). Простейший метод — это действительно метод взвешивания, описанный ivan219.

Кстати, для звукового эквалайзера во многих случаях выгодно использовать биквадратные IIR-фильтры, а не линейнофазовый FIR, получающийся методом взвешивания. Вычислительные затраты обычно получаются значительно меньше.

Если это комбинация устройств "конволюционный процессор + эквалайзер" работающий в реальном времени, то IIR фильтры будут вносить дополнительную нагрузку на CPU хотя их параметры и можно изменять во времени, а FIR- фильтр получится автоматически после предварительного преобразования исходной ИХ - все равно свертку надо делать, поскольку это основная задача. Да, кстати, где можно почитать про другие методы, кроме метода взвешивания?(поскольку мне тоже интересна данная тема)

[Alexey Lukin]

Если это комбинация устройств "конволюционный процессор + эквалайзер" работающий в реальном времени, то IIR фильтры будут вносить дополнительную нагрузку на CPU хотя их параметры и можно изменять во времени, а FIR- фильтр получится автоматически после предварительного преобразования исходной ИХ - все равно свертку надо делать, поскольку это основная задача.

Тут всё зависит от соотношения длин имеющейся свёртки и нового эквалайзера, т.к. добавление эквалайзера удлиняет исходное ядро свёртки. Для руления низких частот приходится использовать FIR-фильтры длиной несколько тысяч или даже десятков тысяч отсчётов. А если делать IIR, то зачастую хватает нескольких умножений на входной отсчёт.

Да, кстати, где можно почитать про другие методы, кроме метода взвешивания?(поскольку мне тоже интересна данная тема)

Я изучал это по Рабинеру/Голду, но там ниже есть даже отдельная книжка: "Digital Filter Design".

[Mad_max]

2 Alexey Lukin Вы не верите, что сигнал можно обрабатывать в частотной области

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

[Serg76]

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

Ну и зачем Вы мне это поясняете, я это и так прекрасно знаю Или Вам просто скучно?

Сообщение отредактировал Serg76 - Feb 5 2011, 11:36

[Mad_max]

Ну и зачем Вы мне это поясняете, я это и так прекрасно знаю Или Вам просто скучно?

да я и не Вам это :-) а тому кто не верит в перемножение спектров

[petrov]

да я и не Вам это :-) а тому кто не верит в перемножение спектров

Проблема как раз в том что часто верят вульгарно, типа взял FFT от блока, занулил ненужные частоты, обратное FFT и вот готов хороший перестраиваемый фильтр.