|
|
Я делаю фильтрацию так.
Введем обозначения.
1. Пусть H - частотная характеристика фильтра длиной m*N (она является желаемой).
2. h = real(ifft(H)) - импульсная характеристика желаемого фильтра длиной m*N.
3. Разбиваем h на m частей. Получаем m частей длиной N. Обозначим их h1, h2, ..., hm.
4. Дополняем каждую из частей нулями, количество которых также N. Обозначим полученные части символами hi0:
hi0=[hi zeros(1, N)]; Эти части имеют длину 2*N.
5. Обозначим Hi0 = fft(hi0); - частотная характеристика длиной 2*N кусочка фильтра.
Все части Hi0 мы рассчитали один раз в матлабе и запихнули их в память ДСП.
Теперь фильтрация.
6. Пусть входной сигнал x, который надо пропустить через фильтр с ИХ h, разбит на блоки длиной N.
Обозначим эти части: x1, x2, x3, ...
7. Будем фильтровать блоками. Т.е. выходной сигнал y мы также должны получить в виде блоков длиной N.
Обозначим их: y1, y2, y3, ...
8. Фильтруем так: берем сигнал x10=[x0 x1] длиной 2*N и выполняем от него БПФ: X10 = fft(x10); x0 - нулевой блок.
9. Y11 = X10.*H1; Y12 = X10.*H2; ... Y1m = X10.*Hm;
10. Y1 = Y11; y1 = real(ifft(Y1)); - готовый отсчет.
11. Повторяем шаг 8 для следущего сигнала: x20 = [x1 x2] длиной 2*N и выполняем от него БПФ: X20 = fft(x20);
12. Y21 = X20.*H1; Y22 = X20.*H2; ... Y2m = X20.*Hm;
13. Y1 = Y12+Y21; y2 = real(ifft(Y2)); - готовый отсчет.
14. Повторяем шаг 8 для блока сигнала k 16. Yk = Y1k+Y2(k-1)+Y3(k-2)+...+Yk1; yk = real(ifft(Yk)); - готовый отсчет. Замечу, что для фильтрации мы используем куски ИХ длиной 2*N - т.е. куски вдвое большей длины, чем блок входного сигнала. Это делается, как известно, для того, чтобы избежать последствий циркулярности свертки, даваемой БПФ. Думаю, очевидно, что такой способ является более эффективным, нежели пропускать блок сигнала длиной N через весь фильтр длиной m*N. Проблема. Как быстро получить H0, имея лишь H? Либо есть иной способ вычисления линейной свертки в частотной области без добавления нулей к ИХ. Я думаю, я не один такой, у кого возникала такая задача. По крайней мере лично я с ней сталкиваюсь уже второй раз. Первый раз решил тупо: ifft, добавление нулей, fft. Т.е. также, как написал выше. Тогда процессорных ресурсов было достаточно. Но сейчас все хуже. У меня не один, а 14 таких фильтров и все они должны иметь возможность ПЛАВНОЙ (!) перестройки. Но про перестройкой утомлять не стану. На крайний случай у меня есть вариантик, как выйти из данной неприятной ситуации, хотя этот вариант и не самый приятный, но он приемлемый. И все это возникает из-за того, что надо к ИХ нули добавлять! Без этого все как по маслу бы прошло.
E-mail:
И так далее.
Но это даже не так важно. Важно то, что в любом случае нам для процедуры фильтрации с использованием fft необходимо иметь частотную характеристику фильтра (или всех его кусков), полученную взятием БПФ от ИХ полученной из ИХ исходного фильтра путем добавления к ней нулей (см. шаги 4 и 5).
Вот здесь проблема и возникает. Она такова. Мы рассчитали фильтр H прямо в частотной области. Но для того, чтобы начать фильтрацию (в той же частотной области), нам не достаточно иметь H. Нам нужно: в простейшем случае, когда m=1, сделать H0 = fft( [real(ifft(H)) zeros(1, N)]; И вот уже используя H0, мы можем делать overlap-add или overlap-save.
Если же разбивать ИХ на куски, то получится более сложный алгоритм, который я привел выше.
Ответы
Перейти к списку ответов
|||
Конференция
|||
Архив
|||
Главная страница
|||
Содержание
info@telesys.ru
