Resampling 44100 -> 12800, нужна помощь по реализации. Опции

[evil_laugh]

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

Прошу помочь мне с алгоритмом, который бы осуществлял изменение частоты дискретизации сигнала с 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

[evil_laugh]

Serg76, похоже, Вы познали Дао. У меня только один вопрос: КАК? Можно, пожалуйста, подробнее про алгоритм?

[Serg76]

Serg76, похоже, Вы познали Дао. У меня только один вопрос: КАК? Можно, пожалуйста, подробнее про алгоритм?

да как и описывал выше на основе полифазника, здесь ничего сверхестественного, правильно подберите нужное окно. я взял распространенное Наттала, правда пришлось все-таки увеличить порядок фильтра до 300, вроде как в -70 дБ влезли (пардон, -60 дБ) А вообще, -60 дБ - это еще по-Божески, где-то попадалось на просторах Инета ресамплер как раз звукового диапазона с подавлением чуть ли не -120 дБ.

Сообщение отредактировал Serg76 - Dec 20 2012, 19:02

[evil_laugh]

да как и описывал выше на основе полифазника, здесь ничего сверхестественного, правильно подберите нужное окно. я взял распространенное Наттала, правда пришлось все-таки увеличить порядок фильтра до 300, вроде как в -70 дБ влезли (пардон, -60 дБ)

То есть Вы делаете простую однокаскадную передискретизацию, но с таким вот хитрым фильтром? А на какую частоту дискретизации фильтр считался? Я так понимаю, Вы не делали "честное" повышение частоты дискретизации до 44100*128, а сделали как-то по-другому? А как?

[Serg76]

То есть Вы делаете простую однокаскадную передискретизацию, но с таким вот хитрым фильтром? А на какую частоту дискретизации фильтр считался? Я так понимаю, Вы не делали "честное" повышение частоты дискретизации до 44100*128, а сделали как-то по-другому? А как?

Да, однокаскадная децимация. В полифазной фильтрации нет ничего "хитрого", она известна давно, есть даже статья в википедии http://ru.wikipedia.org/wiki/%D0%9F%D0%B5%...%86%D0%B8%D1%8F
довольно понятно все описано. никакого повышения/понижения делать не надо. преимущество в том, что все вычисления производятся на выходной частоте дискретизации.

[evil_laugh]

Да, однокаскадная децимация. В полифазной фильтрации нет ничего "хитрого", она известна давно, есть даже статья в википедии http://ru.wikipedia.org/wiki/%D0%9F%D0%B5%...%86%D0%B8%D1%8F
довольно понятно все описано. никакого повышения/понижения делать не надо. преимущество в том, что все вычисления производятся на выходной частоте дискретизации.

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

[Serg76]

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

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

[alex_os]

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

Т.е метод эквивалентный предложенному ув. termit ? А сколько фаз и какая длина каждой фазы ?

[Serg76]

Т.е метод эквивалентный предложенному ув. termit ? А сколько фаз и какая длина каждой фазы ?

Насколько я понимаю, thermit предлагал Лагранжа с последующим использованием в структуре Фарроу. Полифазная фильтрация - это несколько другое. В принципе, можно попробовать еще снизить уровень искажений за счет применения квадратичной интерполяции вместо линейной, но это создаст существенную дополнительную нагрузку на процессор. А так сейчас использую 8 фаз х ~300 коэффициентов каждая, но это с избытком по уровню подавления, так сказать, на скорую руку. При желании число коэффициентов можно уменьшить.

Сообщение отредактировал Serg76 - Dec 21 2012, 07:00

[evil_laugh]

Насколько я понимаю, thermit предлагал Лагранжа с последующим использованием в структуре Фарроу. Полифазная фильтрация - это несколько другое. В принципе, можно попробовать еще снизить уровень искажений за счет применения квадратичной интерполяции вместо линейной, но это создаст существенную дополнительную нагрузку на процессор. А так сейчас использую 8 фаз х ~300 коэффициентов каждая, но это с избытком по уровню подавления, так сказать, на скорую руку. При желании число коэффициентов можно уменьшить.

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

1. Нужно точно знать, какую фазу (точнее, две фазы, или больше, если использовать Лагранжа порядком выше) использовать при вычислении каждого из 128 выходных отсчётов. Либо вычислять это на месте (занимает время), либо заранее рассчитыватт и загонять в lookup table (минимум +256 байт, максимум - зависит от лагранжа).
2. Нужно точно знать, в каком месте между двумя вычисляемыми фазными отсчётами находится интерполируемый, для каждого из 128 интерполируемых. Минимум +256 байт, ну или опять же довольно сложная арифметика на месте.
3. Нужно точно знать, сколько новых отсчётов нужно вдвинуть в буфер при вычислении каждого из 128 отсчётов. Опять либо время, либо +128 байт.
4. Ещё нужно учитывать момент, когда интерполируемый отсчёт находится между последним фазовым и следующим за ним (как бы "нулевым"). Для вычисления "нулевого" фазового нужно временно вдвинуть в буфер ещё один входной отсчёт, профильтровать и снова вернуть буфер на место. (Блин, коряво объяснил, но опытные товарищи, надеюсь, поймут). А это проверка условия на каждом выходном отсчёте. Опять же, время.

Вот так как-то... Всё это работает, но не так эффективно, как хотелось бы К тому же при таком количестве фаз и коэффициентов я не укладывался в требования по подавлению. "По-правильному" фаз должно быть 128, как же оно с восемью так хорошо работает?



Сообщение отредактировал evil_laugh - Dec 22 2012, 08:32