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]

А если так попробовать:
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

[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);

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