Программный фильтр низких частот

stm7543347

Дано: функция f: R -> R, у которой, в принципе, можно подсчитать значение в любой точке (на практике это отображение из компьютерного типа double в него же).
Нужно: получить g, которая как можно более точно совпадает с f на частотах до FA и уходит как можно сильнее в ноль на частотах выше FN, FA < FN. Ну да, дискретизировать хочу.
Чего сделать? Мне пока приходит в голову только веселый бред, навроде: взять много точек слева и справа от данной и по этим точкам устроить какбе локальную дискретную свертку f с sinc. Поскольку вычислительные методы я не разрабатывал никогда, более умного ничего не приходит. :crazy:
ЗЫ. А без фильтра получаются такие вкусные басы от алиазинга в пиле и прямоугольнике. :dj: :grin:

reptilia

а что происходит между F_a и F_n?

stm7543347

Спад, разумеется. Неважно, какой. По идее, F_a - предел слышимости.

reptilia

в качестве гладкого фильтра можно взять функцию ферми-дирака, например.
g(x) = w(x)*f(x x - частоты, w - функция ф-д
но почему бы не использовать
"double g(double x) { return x < f_a ? f(x) : 0.0; }"?

natunchik

Делаешь преобразование Фурье, обнуляешь верхние частоты, делаешь преобразование Фурье.
Если влом: напускаешь на неё фильтр вида r[n] = q * x[n] + (1 - q) * r[n - 1], значение q берешь из http://en.wikipedia.org/wiki/Lowpass_filter

vvasilevskiy

Делаешь преобразование Фурье, обнуляешь верхние частоты, делаешь преобразование Фурье.
++1
Только Фурье быстрое лучше брать

reptilia

:) если f задана во времени - то свертка с функцией ф-д вместо умножения
п.с. пожалуй, это будет быстрее двукратного преобразования фурье

vvasilevskiy

если f задана во времени - то свертка с функцией ф-д вместо умножения
п.с. пожалуй, это будет быстрее двукратного преобразования фурье
Неа, наоборот-свертка-двойной цикл, а бпф-полуторный

reptilia

посмотрел, пишут, что быструю свертку считают через быстрое преобразование фурье.
т.е. и fft, и быстрая свертка считаются примерно одинаково быстро - за O(n log n).

vvasilevskiy

Во, я это и имел в виду а n log n-назвал полуторным циклом

natunchik

А если фурье делать и жестоко обнулять, не получится артефактов разных разве? Плюс, он ведь по дефолту на отрезке работает, если брать за отрезок всю интересующую длительность может больно жирно получиться, а бить на отрезки тоже нетривиальная задача!
Мне всё больше и больше нравится фильтр первого порядка, на самом деле. Глушит частоты экспоненциально сильно, быстрый, тривиальный в реализации.

stm7543347

Делаешь преобразование Фурье, обнуляешь верхние частоты, делаешь преобразование Фурье.
Это то же самое, что свернуть с синком.
Так что, если хотите, могу спросить: как устроить свертку либо преобразование Фурье?

stm7543347

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

stm7543347

Подумал, что вместо синка можно взять фильтр Ланцоша, но, по-моему, сравнительная сладость хрена и редьки на нижеприведенной диаграмме.

stm7543347

Мне всё больше и больше нравится фильтр первого порядка, на самом деле.
Но он не подходит для случая FN / FA = 1,2. Это если рендерить цифровой сигнал на 48 кГц, считая порогом слышимости 20 кГц. Фильтр первого порядка дает 6 дБ на октаву, а при таком соотношении частот ты даже не заметишь, что он был. Алиазинг весь останется.

natunchik

А если его десять раз применить?

stm7543347

Ну да, я и сам об этом подумал.
Потом подумал, что при этом десять раз повернется фаза (фаза при линейной фильтрации поворачивается пропорционально частоте). Уже, возможно, на слух будет заметно.
Однако все это не то. :crazy:
Дело-то вот в чем.
Какие мы берем точки в последовательности? Логично предположить, что те же самые, в которых мы будем брать сэмплы.
Мы ставим на прямой равномерную сетку точек, берем последовательность f(xn проводим фильтрацию и... все это было зря.
Потому что в тот момент когда мы выбрали последовательность и посчитали значения функции в ее узлах... мы уже дискретизировали функцию - разумеется, без фильтра. Именно с алиазингом. Все дальнейшие ухищрения ни к чему не приведут: паразитные басы, появившиеся от него, никаким фильтром НЧ не изгонишь.
Сделаем хитрее. Возьмем какие-нибудь другие точки, большего разрешения, в них проведем фильтрацию и дискретизируем интерполированный сигнал. Что получилось? Правильно. В момент выбора точек фильтрации опять произошел алиазинг. Паразитные басы потом никуда не денутся.
Черт, это смертельная ловушка! :crazy:
Фильтр в данной задаче должен быть - в виде преобразования функций - такого, что его можно записать аналитически, но при этом посчитать на компьютере с любой произвольной точностью. Я ведь на хаскеле все пишу. :crazy:
Медитируя над этими предметами, понял, что цифровой звук - зло, что звуковые карты нас обманывают, что мир изменился и никогда уже не будет таким юным.
Блин, хоть обратно на CSound переходи. :crazy:

romanenkoroman1

получится фильтр 10го порядка :)
а если так вопрос поставить: непрерывное преобразование Фурье - это интеграл. Соответственно и численные методы применять из этого раздела: функцию фильтра-то небось тоже можно аналитически задать
2й способ: проквантовать всё с высокой частотой дискретизации (с такой, чтобы алаясинг уже слабо докучал а потом уже юзать дискретные фильтры

stm7543347

функцию фильтра-то небось тоже можно аналитически задать
sinc x = sin(pi * x) / (pi * x)

Я вот и думаю в сторону численных интегралов.
с высокой частотой дискретизации (с такой, чтобы алаясинг уже слабо докучал)
А такое бывает? :crazy:

romanenkoroman1

не могу быстро найти, но тебе вроде нужно ботать "аналоговые фильтры", "s-преобразование"
насчёт фазы можешь не переживать: если тебе не в реальном времени считать надо, то симметричные фильтры (они используют значения из "будущего") фазу не нарушают. Простое объяснение: фильтранул сигнал в одну сторону, развернул справа налево, фильтранул, обратно развернул - фаза осталась на месте

romanenkoroman1

синк - говно, в нём корень всего алаясинга, найди что-нибудь гладенькое :) размажь ступеньку :)
амплитуда должна убывать с ростом частоты, разве нет? иначе сигнал получится неизмеримым на отрезке

stm7543347

синк - говно, в нём корень всего алаясинга, найди что-нибудь гладенькое :) размажь ступеньку :)
Ну какбе...
Скажем, фильтр Ланцоша внутри окна имеет вид sinc(x) * sinc(x / a и АЧХ его - трапеция.

romanenkoroman1

за своё 3-месячное баловство с DSP я с этим фильтром не сталкивался :( дальше вряд ли что смогу подсказать

natunchik

Что-то я не очень понимаю, что ты делаешь вообще.
Откуда паразитные басы-то берутся? Причём именно басы? Как тебя lowpass filter от басов может избавить? Что за фигня вообще?
Ты не пробовал обычную синусоиду сквозь свою прожраммку прогонять?

romanenkoroman1

ботай aliasing
берём синусоиду y = sin( pi(A+a)t) A-целое, a-маленькое
квантуем: t=2n
получаем: y(n) = sin ( 2pi(A+a)n) = sin(2pi a n) = sin(pi a t) - из писка получили басы

stm7543347

В общем, уже доступно объяснил. Я еще могу пофлудить, что теорема Котельникова (Найквиста, Шеннона и товарищей) гарантирует восстановление исходного сигнала из PCM-оцифровки, только если в исходном сигнале не было гармоник выше F/2, где F - частота оцифровки (в моей формулировке 2 * FN). Эту частоту в англоязычных интернетах называют частотой Найквиста. Скажем, если ты записываешь звук в пресловутом CD-качестве (16/44100/Stereo тебе надо потрудиться, чтобы в нем не было гармоник выше 22050 Гц, по крайней мере так, чтобы амплитуда их была не выше -48 дБ от максимальной, на которую настроен АЦП.
Любая дискретизация есть, разумеется, неинъективное отображение L2 -> R^N. Поэтому неотфильтрованные сигналы склеиваются с некоторыми отфильтрованными. Часто при этом происходит описанный выше алиазинг. Попробуй банально сгенерировать честную цифровую прямоугольную
***   ***   ***   ***

*** *** ***

или пилообразную волну
   *   *   *   *
* * * *
* * * *
* * * *
, сам увидишь, сколько там говна услышишь, какой вкусный, жЫрный получается бас, где не ждали... :dj:

romanenkoroman1

почему ты всё ботаешь цифровые фильтры? ботай аналоговые, это твой случай. В цифровых проблемы бы не было - двусторонний Butterworth нужного порядка

stm7543347

Я ботаю, как программно реализовать аналоговый фильтр. А вот что значит предыдущая фраза, я все никак не могу придумать.

romanenkoroman1

что-то я подумал, что это нереально, если ты интегралы не сможешь аналитически взять. Если ты их вычисляешь по сетке, то все частоты, укладывающиеся внутрь сетки ты всё равно упускаешь. Поэтому всё, что ты можешь сделать, это оценивать скорость сходимости (в L2) и юзать равенство Парсеваля для погрешности в спектральном представлении.
За звуковые карты не переживай: по уму, там аналоговые фильтры реализованы аппаратно (в виде катушек и конденсаторов поэтому сэмплинг проходит правильно

natunchik

А. То есть у тебя как бы муар получается, да? Высокочастотная периодическая компонента интерферирует с, типа, решёткой дискретизации?
Окей, а кто мешает взять в десять или в сто раз более высокую частоту дискретизации и спокойно убить всё лишнее?
Я так понимаю, если у тебя по-прежнему остаются существенные артефакты при увеличении частоты дискретизации в сто раз, то у тебя что-то очень не так с дискретизируемой функцией (ну или глюк в имплементации). Потому что, на секундочку, энергия пропорциональна частоте умноженной на амплитуду и для звука тоже, так что если у тебя действительно такое происходит, это значит, что твоя функция выдаёт практически всю энергию в неслышимой части спектра. Это странно!

stm7543347

энергия пропорциональна частоте умноженной на амплитуду и для звука тоже, так что если у тебя действительно такое происходит, это значит, что твоя функция выдаёт практически всю энергию в неслышимой части спектра. Это странно!
http://en.wikipedia.org/wiki/Square_Wave
[math]\begin{align} x_{\mathrm{square}}(t) & {} = \frac{4}{\pi} \sum_{k=1}^\infty {\sin{\left 2k-1)2\pi ft \right )}\over(2k-1)} \\ & {} = \frac{4}{\pi}\left (\sin(2\pi ft)+{1\over3}\sin(6\pi ft)+{1\over5}\sin(10\pi ft) + \cdots\right ). \end{align}[/math]
Для любой гармоники частота, умноженная на амплитуду, дает единицу. Энергия, как ты ее определил, именно и получается бесконечной. Не вижу ничего странного.

stm7543347

Если ты их вычисляешь по сетке, то все частоты, укладывающиеся внутрь сетки ты всё равно упускаешь.
А что, если просто взять сетку так, чтобы ошибку численного интегрирования сделать не больше погрешности квантования (каковое будет следующим этапом)?

stm7543347

по уму, там аналоговые фильтры реализованы аппаратно (в виде катушек и конденсаторов поэтому сэмплинг проходит правильно
Меня больше успокоила бы фраза "аналоговые фильтры (в виде катушек и конденсаторов) там реализованы по уму". Во всяком случае, я для себя уже решил, что фильтр в ALC650 и фильтр в Lynx Two - это два разных фильтра.

romanenkoroman1

А что, если просто взять сетку так, чтобы ошибку численного интегрирования сделать не больше погрешности квантования (каковое будет следующим этапом)?

если ты можешь это обеспечить - без проблем, но по-моему это в каком-то смысле эквивалентно сначала сэмплингу с высокой частотой, а потом ресамплингу с дискретным фильтром
у непрерывных (периодических) фунций коэффициенты (ряда) Фурье убывают как 1/n^2 (кажется, на 100% не уверен соответственно с непрерывной 1й производной - как 1/n^3
но когда речь идёт о ступеньках, пилах и прочих разрывных хренях - 1/n. поэтому при 16 бит на сэмпл тебе придётся рассматривать ~2^16 гармоник
непериодические функции (с непрерывным спектром) должны в этом же духе себя вести

stm7543347

Сидел на досуге, считал АЧХ. Получилось, что: если взять преобразование F: f(t) -> f(t Ff(t) есть средневзвешенное значений f в точках, взятых в окрестности t с одинаковым шагом и симметричными весами (то же самое, что прямоугольная аппроксимация свертки с синком произвольной четной функцией в прямоугольном окне АЧХ у такого преобразования - функция периодическая, гы-гы-гы.
Оставить комментарий
Имя или ник:
Комментарий: