Low-pass filter

← →
Yr2 (2004-02-02 14:34) [0]

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

← →
Digitman (2004-02-02 14:41) [1]

в общем случае см. алгоритмы быстрого прямого/обратного преобразования Фурье (FFT - Fast Fourier Transforming)

конкретных реализаций этих алгоритмов в сети навалом

← →
Brahman (2004-02-02 14:47) [2]

Для начала, необходимо определить требования к фильтру:
- полоса пропускания с заданной неравномерностью;
- полоса затухания -"-
- "крутизна" скатов АЧХ
- тип реализации IIR/FIR

Простейший фильтр НЧ (первого порядка)
Затухание -6 дБ на декаду

// Инициализация
// tau - период дискретизации
// Ws = 1/T - частота среза по уровню 0.7
// X - входное значение

B := exp(-tau/T);
A := 1 - B;
Yii := 0;
...

//Фильтрация:
Yi := A*X + B*Yii;
Yii := Yi;

← →
Brahman (2004-02-02 14:48) [3]

Digitman © (02.02.04 14:41) [1]
в общем случае см. алгоритмы быстрого прямого/обратного преобразования Фурье (FFT - Fast Fourier Transforming)

Это слишком:)
В зависимости от условий наверняка можно проще.

← →
Digitman (2004-02-02 15:00) [4]

> Brahman © (02.02.04 14:48) [3]

а завтра наверняка его условия изменятся и потребуют более глубокого анализа/синтеза .. так что лучше сразу реализовать FFT , чтобы потом не было проблем с переделкой программы под новые условия

← →
Yr2 (2004-02-02 15:09) [5]

To Brahman
Требования к фильтру:
1. В идеале АЧХ должна быть Г-образной, причём "Г" повернуто в другую сторону.
2. Частота среза должна задаваться в виде параметра функции фильтра, причём её значение 200-400 Гц.
Правильно ли я понимаю, что, чтобы приблизиться к "идеалу", мне нужно поставить последовательно несколько фильтров 1-го порядка? Или формулы расчёта для фильтра N-го порядка другие?
"тип реализации IIR/FIR" - не оговаривается. Можете уточнить, что это меняет?
И последний вопрос: не совсем пока понял роль функции exp(-tau/T); Можно её заменить "мгновенным обрезанием"?

To Digitman
Не хотелось Вам отвечать, но это сделал Brahman. FFT-варианты я уже просмотрел. Это слишком ресурсоемко. Да и не нужно в задаче получать спектр (пусть промежуточный) из N частот.

← →
TUser (2004-02-02 15:11) [6]

Если проще - тогда использвать какую-нибудь матрицу. Ищи на слова "mathematical morphology"

← →
Digitman (2004-02-02 15:22) [7]

> Yr2 (02.02.04 15:09) [5]

> FFT-варианты я уже просмотрел. Это слишком ресурсоемко

как будет угодно.
наше дело - предложить ..

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

к тому же вынужденное каскадирование Фильтров 1-го порядка (для достижения требуемой результирующей АЧХ) будет ощутимо уменьшать производительность алгоритма такого каскадного фильтра в целом, приближая ее к производительности FFT-алгоритма

← →
Brahman (2004-02-02 16:03) [8]

С "идеалом" все понятно, но идеального не бывает, как мы знаем.

Для практических применений вполне бывает достаточно фильтра
4-6 порядка (2-3 фильтра второго порядка, если делать какадную реализацию)

Фильтр первого порядка приведен как простейший случай.
Это RC-цепь с экспонентциальной переходной характеристикой и поэтому там exp(..)

Выбор IIR/FIR - меняет сложность реализации/расчета и еще некоторые параметры

IIR (по старому БИХ-фильтр или фильтр с "бесконечной" памятью)
Для реализации ФНЧ второго порядка:

Result := fGain * (fB0*(fValue + 2*Xi + Xii) - (fA1*Yi + fA2*Yii));

Т.е. обязательных 4 сложения и 4 умножения.
Для фильтра 6 -го порядка (каскадное соединение 3х ФНЧ второго порядка) имеем: по 12 сложений и умножения.
Никакой FFT такого сравнения не выдержит:)

FIR-фильтры (КИХ-фильтры) более ресурсоемки, тк реализуют операцию свертки y = sum(Ai*Xi)
Для создания FIR-фильтра, аналогичного по фильтрующим свойствам IIR-фильтру шестого порядка достаточно i = (100).
Т.е. это будет по 100 сложений/умножений.

Если операция производится не в real-time, то это все мелочи.

← →
Brahman (2004-02-02 16:08) [9]

Инициализация IIR-фильтра:
vFc - частота среза, Гц
vFs - частота дискретизации, Гц
vDzeta - параметр "колебательность" и при 0.707 - стандартный фильтр Баттерворта
vGain - коэф.усиления

Procedure TFilterBQ.SetFilter(vFc,vFs,vDzeta,vGain: double);
var d1,d2,denom: double;
begin
Clear;
fFc := vFc; fFs := vFs; fDzeta := vDzeta; fGain := vGain;
fKs := 0.5/tan(PI*fFc/fFs);
d1 := 4*fDzeta*fKs;
d2 := 4*fKs*fKs;
denom := 1 + d1 + d2;
fB0 := 1/denom;
fA1 := 2*(1 - d2)/denom;
fA2 := (1 - d1 + d2)/denom;
end;

← →
Yr2 (2004-02-02 16:32) [10]

To Brahman:
Спасибо. Я сделал фильтр первого порядка как Вы написали в [2].
Это работает. Единственно, что изменил присвоение Yii:=0; на Yii:=Massiv[0] - это не искажает начало сигнала.
Уже вижу, что одного каскада маловато. Попробую Ваш код [9].
Уточните, пожалуйста, что есть vDzeta? как его оценить и какие из практики ставят значения?

← →
Brahman (2004-02-02 17:02) [11]

dzeta - параметр "колебательность"
При 0.707 (sqrt(2)) - получается фильтр Баттерворта
(максимально плоская полоса пропускания)

← →
Yr2 (2004-02-02 17:58) [12]

Спасибо, Brahman. Ваш алгоритм работает. Но я наблюдаю небольшие странности. Наверное я неправильно понял некоторые переменные.
Если, как Вы пишете, dzeta = sqrt(2), то это не 0.707, а 1.41.
И не могли бы Вы точно написать, как в формуле
Result:=fGain * (fB0*(fValue + 2*Xi + Xii) - (fA1*Yi + fA2*Yii));
вычисляются переменные Yi, Yii, Xi, Xii. Совпадает ли значение Xi с fValue? Или это предыдущий отсчёт?

← →
Brahman (2004-02-02 19:15) [13]

Не углядел: dzeta = 1/sqrt(2) = 0.707...

fValue, X - входные отсчеты, Y - выходные.

fValue == X(i) - текущий
Xi == X(i-1)
Xii == X(i-2)
Result == Y(i) - текущий
Yi == Y(i-1)
Yii == Y(i-2)

← →
Brahman (2004-02-02 20:06) [14]

Если будете каскадировать ФНЧ, то для правильного расчета необходимо учесть следующее
Коэффициент передачи одиночного фильтра на частоте среза Ws составляет 1/sqrt(2) == -3 dB
Если будем использовать трехзвенный фильтр, то на частоте Ws
получим общий коэффициент передачи
Ks = (0.707..)^3 = 0.35
Для того, чтобы обеспечить на частоте Ws общий коэф.передачи
0.707.. (-3дБ) необходимо увеличить частоту среза W1s одиночного
фильтра, так чтобы на частоте Ws каждое звено имело бы коэффициент передачи
(sqrt(2))^(1/3) = 0.8910

АЧХ фильтра Баттерворта
H(W) = 1/(sqrt((W/W1s)^2*n+1))
При n=2(порядок одиночного фильтра) и известных W = Ws и H(W)=0.8910 можно найти частоту среза одиночного фильтра W1s = sqrt(2)*Ws = 1.414*Ws

← →
Yr2 (2004-02-02 20:43) [15]

To Brahman
Всё, большое спасибо! Вы меня удовлетворили :) Точнее, Ваш фильтр 2-го порядка :)
Теперь отфильтрованный сигнал выглядит так как и должен выглядеть.

P.S. Жаль что в форуме нельзя бутылку пива передать...

Low-pass filter Найти похожие ветки