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

[Serg76]

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

[evil_laugh]

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

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

[beaRTS]

Прошу помочь ... бла бла бла

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

[evil_laugh]

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

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

[Serg76]

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

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

[des00]

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

а) пропускание частоты сигнала в полосе от 0 до 6 кГц с неравномерностью АЧХ не более 3 дБ
б) подавление сигнала в полосе 6200 Гц - 6400 Гц не меньше, чем на 60 дБ

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

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

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

[Serg76]

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

[evil_laugh]

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

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

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

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

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

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

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

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

[Serg76]

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

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

[des00]

А Фарроу не будет искажения в сигнал вносить засчёт алиасинга?

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

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

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

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

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

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

Due to the fact that this filter belongs to the set of wave digital filters, it is stable for any coefficients. Moreover, because of its simple and single coefficient its complexity is very small. Such an coefficient can be implemented using a single adder which adds the multiplicand shifted to 2 and 3 bits right.

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

Эскизы прикрепленных изображений

 

[Самурай]

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

Как еще один вариант ("честный") - каскад дробных полифазных КИХ фильтров, например в такой последовательности : 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

[Serg76]

Фильтры в каскадах 4/7 и 4/3 нужны совсем небольшие, 12-16 порядка навскидку.

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

[evil_laugh]

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

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

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

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

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

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

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

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

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

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

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

[Самурай]

Соврал немного про каскады 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 порядка.

[Serg76]

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

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

ну потому как порядок фильтра и будет определяться отношением частот 44100/12800, т.е. оверсемплингом + подбираем нужное окно. А Вы пытаетесь сначала провести целочисленную интерполяцию в 128 раз. конечно при таком оверсемплинге и к-ту подавления -60 дБ у Вас получится нереальный порядок фильтра. АЧХ смотрел на "глаз", поэтому могут быть небольшие отклонения в количестве коэффициентов.

Соврал немного про каскады 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 порядка.

я на это как раз и обратил внимание, но проверять не хочется

[alex_os]

А если так попробовать:
1. Делаем ФHЧ тактовой 44100 с подавлением -60 дБ на 6.2 кГц
2. Полиноминальный интерполятор 3го порядка, или вообще линейный.

Порядок ФНЧ у меня получается что-то около 350. Можно его реализовать в частотной области на основе FFT512,
либо способом частотной маски.

[evil_laugh]

А если так попробовать:
1. Делаем ФHЧ тактовой 44100 с подавлением -60 дБ на 6.2 кГц
2. Полиноминальный интерполятор 3го порядка, или вообще линейный.

Порядок ФНЧ у меня получается что-то около 350. Можно его реализовать в частотной области на основе FFT512,
либо способом частотной маски.

Вот тут не совсем понял... Вы предлагаете сперва ограничить по частоте исходный сигнал, а потом куда-то его интерполировать?..

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

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

Сделал я Фарроу. Кубический. Симметричный. Децимирует с 44100 до 38400. Тестировал так: генерировал испытательный синус 0 dBFS с ЛЧМ модуляцией (от 0 до 20 кГц в течение 100 секунд) и прогонял через наскоро набросанную на Си программку. Результат печален. На частотах до 2-3 кГц всё неплохо, гармоники в сигнал не лезут. С повышением частоты начинаются весёлости. Какие - красочно иллюстрирует рисунок ниже. На входе синус 6 кГц, Fд = 44100. На выходе - вот что.



Хорошо видно, что из-за полиномиальной интерполяции в сигнал лезут ложные гармоники просто из-за ошибки интерполяции. Сгенерировал в маткаде идеальный дискретный синус и нашёл разность полученного мной сигнала и идеального. Полная шкала от -32768 до +32767.



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

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

[thermit]

Результат соответствует действительности.

точность лагранжа повышается повышением частоты дискретизации.
например:

44100->132300 фнч с заданной х-кой (~1300x16 коэффициентов + ~400x16 память фильтра) -> 12800 при помощи лагранжа.
Если 1-го порядка (линейная интерполяция) подавление > =60дб ~10-11 mips.

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

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

[alex_os]

Вот тут не совсем понял... Вы предлагаете сперва ограничить по частоте исходный сигнал, а потом куда-то его интерполировать?..

Понижать частоту дискретизации до требуемой. Тем же фарроу или линейным интерполятором.
Суть в том, что нужно сначала отфильтровать сигнал до требуемых 6.2 кГц а потом переходить
на другую частоту дискретизации. В противном случае все что выше 6.2 кГц неизвестно куда завернется.

upd: По моим оценкам полинома 3го порядка должно хватить чтобы уйти с 44100 на 12800.

Код

Fs = 44100;
K =17   % коэф. передискритезации вверх не важно какой
h0 = ones(1,K);
h1 = conv(h0, h0);  % интерполяция полиномом 1го порядка
h1 = h1/sum(h1);
h2 = conv(h1, h1);
h2 = h2/sum(h2);   % интерполяция полиномом 3го порядка
hfir = firpm(380, [0,6000/(Fs/2), 6200/(Fs/2),1], [1,1,0, 0], [0.01,1]);
figure(1);
freqz(h2,1);
F = [1:Fs/2];
H = freqz(hfir, 1, F, Fs);
figure(1);
plot(F, 20*log10(abs(H)));

htmp = zeros(1, length(hfir)*K);
htmp(1:K:end) = hfir;

hall = conv(htmp, h2);

figure(2);
freqz(hall,1);

Эскизы прикрепленных изображений

 

[des00]

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

сделал модель в симулинке с интерполятором 4 го порядка, не очень получается, особенно выше 3КГц. На выходном синусе паразитная АМ, возникающая из-за переменной фазы при децимации с 51200/38400 на 12800. Что-то мне подсказывает что на полиномиальных интерполяторах лучше не сделать, ну разве что интерполировать сначала раз в 16, а потом опускаться.

[alex_os]

сделал модель в симулинке с интерполятором 4 го порядка, не очень получается, особенно выше 3КГц. На выходном синусе паразитная АМ, возникающая из-за переменной фазы при децимации с 51200/38400 на 12800. Что-то мне подсказывает что на полиномиальных интерполяторах лучше не сделать, ну разве что интерполировать сначала раз в 16, а потом опускаться.

Что-то не так. Берем синус 3 кГц Fs = 44100 и интерполятором 4го порядка напрямую опускаем на Fs=12800.
Должно идеально получаться...

[thermit]

alex_os:
По моим оценкам полинома 3го порядка должно хватить чтобы уйти с 44100 на 12800.

Для 3-степени погрешность интерполяции синуса единичной амплитуды:

err<= 3/8*(pi*f0/fd)^4

fd=44100
для f0=3000 err<= 7.822785453414741e-004
для f0=4000 err<=0.002472386513672
для f0=6000 err<= 0.012516456725464
Топикстартеру это не понравится.

Если выполнить хотя бы 2-х кратный оверсамплинг 44100 -> 88200
для f0=3000 err<= 4.889240908384213e-005
для f0=4000 err<=1.545241571044886e-004
для f0=6000 err<= 7.822785453414741e-004

Это уже топикстартера должно устроить. Порядок фильтра-интерполятора ~700.

[evil_laugh]

Большое всем спасибо за грамотные и содержательные ответы. Картина значительно прояснилась. Попробую повысить частоту дискретизации и потом использовать Лагранжа.

Надо учиться пользоваться Матлабом, а то симуляция различных вариантов много времени занимает

[Serg76]

Сделал я Фарроу. Кубический. Симметричный. Децимирует с 44100 до 38400. Тестировал так: генерировал испытательный синус 0 dBFS с ЛЧМ модуляцией (от 0 до 20 кГц в течение 100 секунд) и прогонял через наскоро набросанную на Си программку. Результат печален. На частотах до 2-3 кГц всё неплохо, гармоники в сигнал не лезут. С повышением частоты начинаются весёлости. Какие - красочно иллюстрирует рисунок ниже. На входе синус 6 кГц, Fд = 44100. На выходе - вот что.

Накидал и я свой интерполятор по алгоритму, который предлагал выше. Децимирует с 44100 до 12800. Тестируемый сигнал - бигармонический сигнал с частотами 3500 и 14000 Гц. Вот что получилось:

до



после



вроде как требования выполняются

[alex_os]

Для 3-степени погрешность интерполяции синуса единичной амплитуды:

err<= 3/8*(pi*f0/fd)^4
...

Похоже на то. Я сжульничал и вместо лагранжа 3ей степени вставил нечто другое. Правильные рога для Лагранжа
1, 3, 5 порядков такие (см. в самом низу). Выходит что нужен Лагранж 5 степени чтобы с 44100 перейти на 12800 при <-60 дБ наложений.

Код

Fs = 44100;
K =11   % коэф. передискритезации вверх не важно какой

hfir = firpm(380, [0,6000/(Fs/2), 6200/(Fs/2),1], [1,1,0, 0], [0.01,1]);
F = [1:Fs/2];
H = freqz(hfir, 1, F, Fs);
figure(1);
plot(F, 20*log10(abs(H)));

for p=1:2:5
    h = intfilt(K,p,'Lagrange');
    h = h/sum(h);
    
    htmp = zeros(1, length(hfir)*K);
    htmp(1:K:end) = hfir;

    hall = conv(htmp, h);

    figure(2);
    hold on;
    freqz(hall,1);
end

Накидал и я свой интерполятор по алгоритму, который предлагал выше. Децимирует с 44100 до 12800. Тестируемый сигнал - бигармонический сигнал с частотами 3500 и 14000 Гц. Вот что получилось:

Частоты неинтересные у Вас . Интересно что получится если подать два синуса 6 кГц и 6.2 кГц...

Эскизы прикрепленных изображений

 

[Serg76]

Частоты неинтересные у Вас . Интересно что получится если подать два синуса 6 кГц и 6.2 кГц...

пожалуйста

[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, а затем для каждого элементарного фильтра пересчитываем эту характеристику с учетом смещения. в процессе фильтрации для каждого выходного отсчета вычисляем отклики двух соседних фильтров и проводим линейную интерполяцию между ними. поэтому выходной спектр и "обогащается" гармониками, которых не было во входном сигнале.

[evil_laugh]

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

Вот теперь мозаика сошлась в единое целое. Спасибо. Я пробовал делать подобное (пункт 4 первого сообщения темы), правда, с другой стороны подошёл немного. Попробую сделать и так, как Вы сказали. Заодно время выполнения алгоритма и общую требуемую память подсчитаю.

[alex_os]

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

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

[Serg76]

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

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

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

[thermit]

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

1. Фарроу - один из способов организации вычислений интерполяции лагранжа. Достоинство - минимум умножителей. Поэтому просто лагранж.

2. Вы тоже используете интерполяцию лагранжа. 1-го порядка. Впрочем, таковой она будет для любой интерполяции 1-го порядка. Интерполяция к-тов фильтра полностью эквивалентна интерполяции непосредственно сигнала.

3. Вы используете 8-кратный оверсамплинг (8 фаз фильтр 2399-го порядка)

4. Можно фильтр покороче (оверсамплинг поменьше) - но интерполятор порядка повыше.

5. Так что, это именно то, что я предлагал. Что-то принципально новое изобрести здесь крайне сложно. (Вероятно - невозможно).

[Serg76]

1. Фарроу - один из способов организации вычислений интерполяции лагранжа. Достоинство - минимум умножителей. Поэтому просто лагранж.

2. Вы тоже используете интерполяцию лагранжа. 1-го порядка. Впрочем, таковой она будет для любой интерполяции 1-го порядка. Интерполяция к-тов фильтра полностью эквивалентна интерполяции непосредственно сигнала.

3. Вы используете 8-кратный оверсамплинг (8 фаз фильтр 2399-го порядка)

4. Можно фильтр покороче (оверсамплинг поменьше) - но интерполятор порядка повыше.

5. Так что, это именно то, что я предлагал. Что-то принципально новое изобрести здесь крайне сложно. (Вероятно - невозможно).

Интерполяция одна и та же - полиномиальная, не спорю. Схемы реализации разные - либо в виде структуры Фарроу, либо полифазная. Уменьшить количество коэффициентов тоже можно - оптимизация по Ремезу, подбор окна, компромисс между порядком полинома и количеством фаз и т.д.

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

[thermit]

Serg76:
Схемы реализации разные - либо в виде структуры Фарроу, либо полифазная.

Одно другому не мешает. Полифазный фильтр выполняет либо вычисление осчетов соответствующих целым фазам, тогда промежуточный отcчет вычисляется посредством интерполяции по вычисленным целофазным отсчетам, либо вычисление непосредственно промежуточного отсчета, но с импульсной х-кой интерполированной по к-там. В обоих случаях нужны:
a) Полиномиальная интерполяция
б) Фильтрация
Теоретически, оба способа эквивалентны. А вот с точки зрения ресурсов - нет.
Пример:
Число к-ов в фазе N линейная интерполяция
Интерполяция по отсчетам:
2 отсчета - 2*N macc
интерполяция 2 macc
итого 2*N+2 macc

Интерполяция к-ов:
интерполяция к-ов 2*N macc
Вычисление отсчета N macc
итого 3*N macc

Так, что способ 1, очевидно, эффективнее.

[Serg76]

Одно другому не мешает. Полифазный фильтр выполняет либо вычисление осчетов соответствующих целым фазам, тогда промежуточный отcчет вычисляется посредством интерполяции по вычисленным целофазным отсчетам, либо вычисление непосредственно промежуточного отсчета, но с импульсной х-кой интерполированной по к-там. В обоих случаях нужны:

Да, да, да и еще раз ДА, по сути одна схема трансформируется в другую. По поводу эффективности: мне, кажется, что это тема более глубокого исследования, просто я привык пользоваться полифазниками, Лагранжем + Фарроу пользовался единожды, тогда мне это не понравилось, может стоит пересмотреть свои взгляды

[Serg76]

Sorry, что-то глюкнуло

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

[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

[thermit]

evil_laugh:
Тут столкнулся с рядом проблем, которые тоже увеличивают время выполнения и требуемую память.
...

Странные проблемы.
Да. Текущую фазу надо вычислять.
По целой ее части выбирается фаза полифазного фильтра. Дробная часть определяет интерполяцию.
Памяти нужно N под коэффициенты фильтра и N/OVERSAMPLING под память фильтра. ~ 2*N/OVERSAMPLING macc на выходной отсчет.
Все вполне эффективно.

ps
С вдвиганием-выдвиганием отсчетов вы перемудрили

[Serg76]

ну давайте, попробуем по пунктам:

1. Нужно точно знать, какую фазу (точнее, две фазы, или больше, если использовать Лагранжа порядком выше) использовать при вычислении каждого из 128 выходных отсчётов. Либо вычислять это на месте (занимает время), либо заранее рассчитыватт и загонять в lookup table (минимум +256 байт, максимум - зависит от лагранжа).

в общем случае необходим параметр, некий аргумент, который будет задавать все моменты, в которых нужно производить соответствующие расчеты, зависит он от оверсемплинга. он определяет момент вычисления выходного отсчета, определяет смещение между двумя интерполируемыми отсчетами и, соответственно, определяет фазы, между которыми находится текущий отсчет. все это вычисляется на лету, т.е. на месте, никаких заранее просчитанных ЛУТов не используется, за исключением таблицы с ИХ фильтров.

2. Нужно точно знать, в каком месте между двумя вычисляемыми фазными отсчётами находится интерполируемый, для каждого из 128 интерполируемых. Минимум +256 байт, ну или опять же довольно сложная арифметика на месте.

определяется аргументом из п.1.

3. Нужно точно знать, сколько новых отсчётов нужно вдвинуть в буфер при вычислении каждого из 128 отсчётов. Опять либо время, либо +128 байт.

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

4. Ещё нужно учитывать момент, когда интерполируемый отсчёт находится между последним фазовым и следующим за ним (как бы "нулевым"). Для вычисления "нулевого" фазового нужно временно вдвинуть в буфер ещё один входной отсчёт, профильтровать и снова вернуть буфер на место. (Блин, коряво объяснил, но опытные товарищи, надеюсь, поймут). А это проверка условия на каждом выходном отсчёте. Опять же, время.

Да есть такой нюанс, Вы правильно заметили. но как-то его надо учитывать.

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

Работает это все очень даже эффективно, если еще оптимизировать количество необходимых коэффициентов, то выглядит это все очень неплохо. кстати очень много ресамплеров звукового диапазона работают именно на основе полифазной фильтрации, какой у них алгоритм используется - я не знаю. вся эта обвязка, отвечающая за синхронизацию различных моментов в ресамплере, очень проста и работает в реальном времени. сам интерполятор (без функций предварительного расчета ИХ фильтров) у меня занимает не более 100 строк кода, ну после соответствующих оптимизаций . При 8-ми фазах это работает, результаты я приводил, можно конечно увеличить число фаз, тем самым снизить погрешность интерполяции, но это приведет к пропорциональному увеличению количества памяти под коэффициенты и незначительному выигрышу С/Ш, поэтому в этом случае я бы лучше перешел от линейной интерполяции к квадратичной, но при этом уже увеличится вычислительная нагрузка, т.к. вычислять уже придется 3 точки вместо 2-х. ну как-то так.