Приветствую уважаемых гуру цифровой обработки

Прошу помочь мне с алгоритмом, который бы осуществлял изменение частоты дискретизации сигнала с 44100 до 12800 Гц. Сигнал - 16битные отсчёты с микрофона, обычный PCM.

Требования к передискретизатору следующие:

а) пропускание частоты сигнала в полосе от 0 до 6 кГц с неравномерностью АЧХ не более 3 дБ
б) подавление сигнала в полосе 6200 Гц - 6400 Гц не меньше, чем на 60 дБ
в) и самое главное - крайне жёсткие требования к объёму используемой оперативной памяти и процессорного времени. Процессор, для которого пишется алгоритм, помимо передискретизации загружен ещё тонной задач. Оперативной памяти для алгоритма может быть выделено не более 2 кБ. Время - как всегда, чем меньше, тем лучше.

Расскажу о том, какие исследования я провёл в этой области.

1. Естественно, что чтобы сделать изменение Fд 44100->12800, нужно сперва повысить частоту в 128 раз, потом понизить в 441 раз. Фильтр, который нужен для такой передискретизации в один каскад, имеет больше 7000 коэффициентов (рассматриваем КИХ-фильтрацию), поэтому вариант с таким алгоритмом отпал сразу.
   
2. Следующим моим шагом было создание многокаскадной системы: сперва повышение частоты в 2*4*16=128 раз, потом её уменьшение в 21*7*3=441 раз. Естественно, использовалось полифазное разбиение фильтров. На входе алгоритма был 441 отсчёт, на выходе - 128. Всё работало и обеспечивало требуемые параметры АЧХ, однако требовалось много места (нужно же где-то хранить 441*2*8 промежуточных отсчётов сигнала на высокой частоте) и много процессорного времени. Поэтому от этого алгоритма тоже пришлось отказаться. Было поставлено требование: только однокаскадная передискретизация.
   
3. Попутно, кстати, была прочитана горка литературы по фильтрам Фарроу, Ланцоша и прочим полиномиальным фильтрам, но ввиду почти полного отсутствия вменяемых реализаций с графиками и прочими данными - я понял только общий принцип и то, что обмануть никого не удастся и всё равно придётся реализовывать огроменный фильтр.
   
4. Дальше попробовал сделать так: разработать фильтр с меньшим количеством отсчётов, а при фильтрации проводить интерполяцию полученных значений. При этом было решено поступиться требованиями к переходной полосе фильтра и сделать так, чтобы на 6200 Гц обеспечивалось подавление не ниже 15 дБ, но все дальнейшие частоты выше какой-то величины (у меня получилось 7 кГц) всё равно должны были подавляться не меньше, чем на 60 дБ. Алгоритм был написан на Си и показал свою работоспособность. Казалось бы, всё в порядке и можно писать программу для DSP, но меня не покидало ощущение того, что можно попробовать что-то ещё. Поэтому плавно переезжаем к...
   
5. ...фильтрации в частотной области, которую я решил рассмотреть как вариант реализации передискретизатора. Что нужно? Правильно, взять ДПФ исходных 441 точек сигнала, занулить ненужные гармоники, выкинуть ненужные 441-128 отсчётов из середины спектра и золотой ключик у меня в кармане. Однако снова возник ряд проблем. Первая - это "как быстро взять ДПФ 441 точек?". Дополнение нулями до 512 точек, судя по всему, не прокатывает - для передискретизации в частотной области придётся выкинуть нецелое количество точек. Ладно, для этого вопроса решение было найдено путём штудирования книжки Блейхута "Быстрые алгоритмы цифровой обработки", откуда я выудил алгоритмы Кули-Таки и алгоритм Гуда-Томаса (который прекрасно работает для числа 441, ибо это 3*3*7*7 - произведение простых чисел). Решение спорное, впрочем, ибо архитектура DSP, ИМХО, ну ни разу не приспособлена для подобного рода выкрутасов со всеми этими адресными тасованиями, тем более в моём дилетантском исполнении. Впрочем, возникает вторая проблема, которая затмевает собой первую: как избежать эффекта Гиббса в данной ситуации? Напрашивается всем известный overlap-save, однако не получится ли так, что время, затрачиваемое на вычисление "лишних" точек, которые потом выкидываются, сведёт весь выигрыш от используемого алгоритма на нет?

Собственно, у меня вопрос к уважаемому сообществу: что можно ещё придумать, чтобы выполнить задачу? Приветствуются пинки ногами, крики "пшёл вон, нуб", а также ссылки на авторитетные источники, личный опыт и другую информацию. Напомню задачу: 44100 -> 12800. пропускание 0-6 кГц с неравномерностью 3 дБ, подавление 6200-6400 минимум 60 дБ, минимум процессорного времени, 2 кБ памяти.

Спасибо


Сообщение отредактировал evil_laugh - Dec 18 2012, 16:03

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

Спасибо за ответ. Фильтры я реализовывал с помощью известного метода замены Ремеза, они вроде как оптимальные по количеству коэффициентов и подавлению/пропусканию. Окна пробовал - результат печален, либо тонна коэффициентов, либо ЧХ от желаемых далеки, как Луна от Земли Может, есть ещё что-то, что я упустил из виду?

если не понятно о чем говорит Serq76 - можно подглядеть в книгу Айфичер Э., Джервис Б. Цифровая обработка сигналов Вроде б, глава 9. ну найдете

Спасибо за ответ. Эту книжку я до дыр зачитал уже Там совсем мало про передискретизацию с нецелым коэффициентом. Моей настольной книжкой в этом отношении стала зарубежная книга Chrochiere&Rabiner "Multirate Digital Signal Processing", но все "лежащие на поверхности" методы оттуда я уже, кажется, испробовал... Хотя не знаю. Вот для того и перечислил всё то, что я уже сделал или изучил, чтобы было понятнее

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

судя по всему в вас оверсамплинг в 44,1/6 = 7,35 раза. Берем кубический фарроу, делаем ресамплинг 44100 -> 38400 (тут только надо посмотреть, куда избыточный шум завернется) ну и затем классический дециматор на 3. Матлаб рисует фильтр 155 порядка (затухание на 6400 60дб, в полосе 6000 неравномерность 3дБ), что для децимирующего фильтра не так уж и много.

ЗЫ. потребуется память под 155 коэффициентов, циклический буфер под 155+3 отсчета данных, MAC и считать каждый третий отсчет %)

UPD. Можно подняться с 44100 до 51200, тогда нужно посчитать куда зеркало встанет (его давить проще ИМХО), ну а дальше децимация на 2 тем же циком или bi-reciprocal фильтрами и фильтрация с возможной коррекцией %)

ну вот, судя по расчетам des00 я вроде не так и далек оказался по части количества коэффициентов.

Вот как раз в уровне подавления собака и порылась. Если уменьшить требования к моему подавлению хотя бы на 10 дБ, то 7 тысяч коэффициентов одним махом превращаются в 3-4. Это явление (правда в "профиль" - уменьшал не подавление, а длительность переходной полосы) было использовано в пункте 4 сабжа. Если делать "по правилам" с соблюдением всех условий - будет так, как есть

А вот такая комбинация мне не приходила в голову. Именно полиномиальный дециматор в каскаде с обычным... Мне нравится идея. А Фарроу не будет искажения в сигнал вносить засчёт алиасинга? Вообще я пробовал записывать речевой сигнал с микрофона, там редко когда за 16-17 кГц частота вылезает, но всплески порой бывают. Они не будут сильно влиять? И что делать, если будут?


Страшные слова... Я таких не слышал даже... Впрочем, с апсемплингом тоже вариант...

Сообщение отредактировал evil_laugh - Dec 18 2012, 16:50

не знаю как Вы там считаете, но только что тоже запустил матлаб и 150 коэффициентов достаточно для обеспечения подавления в 60 дБ.

тут надо смотреть что вы будете делать с помощью Фарроу: интерполировать или децимировать (разница будет в том, как сформировать сигнал дробной задержки). При интерполяции (подняться до 51200) появившиеся зеркала отфильтрует ваш же собственный фильтр на 6КГц. При децимации будет заворот шумов в полосу, но тут надо посчитать куда именно будет заворот, если он попадет вне вашей рабочей полосы, то и пусть будет %) (ну это если у вас там только шум, без каких нить предустановленных узкополосных помех)


в приложении ЧХ интерполятора третьего порядка. На ней виден линейный участок на котором желательно работать. При вашем оверсамплинге рабочей полосе ИМХО ничего не угрожает.

На подобных конструкциях (интерполирующие и децимирующие дробные ресамплеры) маленькую таксу съел, если что обращайтесь помогу моделями в симулинке %)



это простой БИХ фильтр хорошо подходящий для интерполяции и децимации на 2. Его достоинство :

правда подавления в 40дб обычно мало, приходиться ставить 2 звена

Как еще один вариант ("честный") - каскад дробных полифазных КИХ фильтров, например в такой последовательности : 44100*(4/7)*(4/7)*(4/3)*(2/3) = 12800. Основная фильтрация с переходной полосой 6000-6400 выполняется фильтром последнего каскада, это эквивалентно вышеописанному дециматору с 38400 до 12800 и порядком фильтра 155. Фильтры в каскадах 4/7 и 4/3 нужны совсем небольшие, 12-16 порядка навскидку.

Сообщение отредактировал Самурай - Dec 18 2012, 17:10

алиаса не будет при такой "громоздкой" схеме? да и незачем, наверное, делать такой многокаскадный ресамплер при таком низком оверсемплинге. можно напрямую или как советует des00. помнится в одной из задач надо было передисретизировать сигнал в 3000 раз, обошелся 3-мя каскадами .

В таком случае я просто не могу не задать Вам ужасающий своей тупостью и отдающий крайней степенью дилентантизма вопрос: как считаете такой фильтр Вы, что у Вас получается 150 коэффициентов?

Фильтр, который стоит между интерполятором и дециматором при однокаскадном дробном ресамплинге, должен работать на частоте дискретизации, которая получается после интерполяции. В моём случае она равна 44100*128 = 5.625 МГц. Если рассчитать на эту частоту ФНЧ с требуемыми параметрами - получится больше 7 тысяч коэффициентов. Что я делаю не так?


Так-так, чрезвычайно интересно... Завтра попробую провести расчёты... Входной сигнал прост, как три рубля - передаваемая через микрофон человеческая речь. У микрофона АЧХ вряд ли все 20 кГц имеет... Завтра измерю. Спасибо!



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

Спасибо большое. С Вашей помощью собаку, если повезёт, съем и я



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

Соврал немного про каскады 4/7 и 4/3, фильтры там нужны все-таки побольше
Если снова чего не напутал, то получается так:

44100*(4/3) = 58800 - фильтр 12 порядка (подавление 80дБ),
58800*(4/7) = 33600 - фильтр 33 порядка (подавление 80дБ),
33600*(4/7) = 19200 - фильтр 54 порядка (подавление 80дБ)

И последний каскад с основной фильтрацией (про который уже говорилось):
19200*(2/3) = 12800 - фильтр 155 порядка.