2.3. Проектирование полосовых фильтров
Полосовые фильтры данной системы должны обладать хорошими характеристиками фильтрации. После аппроксимации характеристики получившегося цифрового ПФ должны быть максимально приближены к характеристикам идеального ПФ. Фильтр должен иметь резкий спад в переходной полосе, гладкую граничную полосу пропускания и задерживания. Такая точность необходима исходной вычислительной системе для максимального уменьшения вредных шумов и четкого разграничения суммарного сигнала РГМ по полосам.
Выше изложенными характеристиками, на мой взгляд, обладают несколько различных видов фильтров.
Рис.12. Модели АЧХ различных фильтров с порядком n=5
В первую очередь это фильтр Баттерворта, который обладает гладкой переходной полосой пропускания и задерживания, и на высоких порядках имеет резких спад. Существенным плюсом данного фильтра является отсутствие изменения формы АЧХ с увеличением порядка фильтра.
Во вторых это фильтры Чебышева, 1го и 2го типа, обладающих более резким спадом, однако с небольшим биением в полосах пропускания и задерживания соответственно.
И в третьих это фильтр Кауэра (эллиптический), обладающий самым резким спадом из всех перечисленных фильтров. Большим минусом данного фильтра является изменение формы биений и характеристик спада с увеличением порядка фильтра. А также наличие самих биений и в полосе пропускания, и в полосе задерживания.
При проектировании фильтра необходимо учитывать, что требуемая система не должна иметь высокой задержки сигнала, поэтому нецелесообразно использовать фильтры выше 6 или 7го порядка, чтобы не вносить больших дополнительных задержек в блок фильтрации и обработки исходного усиленного сигнала головного мозга.
Для того чтобы РГМ не накладывались друг на друга, нам необходимо четко их разграничивать, поэтому самым рациональным в данном случае будет фильтр Баттерворта 6го или 7го порядка, с четким спадом и отсутствием биений.
Главной частью в проектирование фильтра является создание его передаточной функции. По сути, передаточная функция и есть фильтр. При подстановке дискретных (точечных) значений в которую исходная поступающая величина меняется в зависимости от характеристик фильтра.
Воспользуемся следующим отрезком кода, созданным используя программный пакет MatLab фирмы The MathWorks, Inc. для нахождения передаточных функций с помощью перехода от НФНЧ к ПФ:
Fs1=; % нижняя граничная частота полосы задерживания (Гц)
Fs2=; % верхняя граничная частота полосы задерживания (Гц)
Fp1=; % нижняя граничная частота полосы пропускания (Гц)
Fp2=; % верхняя граничная частота полосы пропускания (Гц)
Ft=1000; % частота отсчетов (Гц)
Rs=19; % Дб
Rp=14; % Дб
ws1=2*pi*Fs1/Ft; % рад/с (нормированные частоты)
ws2=2*pi*Fs2/Ft; % рад/с
wp1=2*pi*Fp1/Ft; % рад/с
wp2=2*pi*Fp2/Ft; % рад/с
Ws1 = 2*Ft*tan(ws1/2); % (граничные частоты аналогового фильтра)
Ws2 = 2*Ft*tan(ws2/2);
Wp1 = 2*Ft*tan(wp1/2);
Wp2 = 2*Ft*tan(wp2/2);
W2p2=Ws1*Ws2/Wp1;
Bw=W2p2-Wp1; %Ширина полосы пропускания
W0=sqrt(W2p2*Wp1); %Центральная частота полосы пропускания
Wst=(Ws2-Ws1)/(Wp2-Wp1); %граничной частоты полосы задерживания
[n, Wn] = buttord(1, Wst, Rp, Rs, 's') %Нахождение порядка НФНЧ
[num, den] = butter(n, Wn, 's');
[BT, AT] = lp2bp(num, den, W0, Bw) % нахождение коэффициентов передаточной функции фильтра
Подставляя необходимые частоты в отрезок кода приведенный выше для полосы пропускания и задерживания, мы можем получить передаточные функции для фильтров, обеспечивающих разделение исходного сигнала по РГМ.
Характеристики фильтра для альфа-ритма
Fs1=8; % нижняя граничная частота полосы задерживания (Гц)
Fs2=14; % верхняя граничная частота полосы задерживания (Гц)
Fp1=8.5; % нижняя граничная частота полосы пропускания (Гц)
Fp2=13; % верхняя граничная частота полосы пропускания (Гц)
Числитель передаточной функции:
103 *(3.81s3 -1.7 * 10-14s2 -1.01*10-11s1-2.2*10-14)
Знаменатель передаточной функции:
1010 *(s6 10-9 +s5 3,12*10-9+ s4 1,33*10-6 + s3 2,7*10-5 + s2 5,6*10-3 +s 5,7*10-2 +7,77)
Порядок фильтра: N = 6
А для построения АЧХ и ЛАЧХ фильтра воспользуемся следующим отрезком кода: (2)
[H, f] = freqz(numd, dend, 256, Ft);
figure(2)
subplot(211)
plot(f, abs(H)), grid
title('Magnitude response')
xlabel('Frequency, Hz ')
subplot(212)
plot(f, 20*log10(abs(H))), grid
title('Logarifmic magnitude response')
axis([0, 11000, -320, 10])
xlabel('Frequency, Hz ')
Рис.13. АЧХ И ЛАЧХ получившегося фильтра для альфа-ритма
Код в Matlab для билинейного преобразования:
[numd, dend] = bilinear(BT, AT, Ft)
Чтобы построить данную АЧХ нам пришлось перейти от аналогового прототипа фильтра к его цифровому варианту, с помощью билинейного преобразования.
Метод билинейного
преобразования основан на использовании
конформного преобразования р-плоскости
с помощью функции
(1)
где
–
положительная константа.
При переходе от ПФ аналогового фильтра к ПФ цифрового фильтра используется следующее выражение:
Соответствие между
частотами аналогового фильтра
и
цифрового фильтра следует из подстановки
в выражение (1)
и,
где
–
интервал дискретизации. При такой
подстановке:
Эквивалентное выражение
Отсюда
Преобразование
является нелинейным, но оно близко к
линейному при небольших .
Частота среза
аналогового прототипа соответствует
частоте среза
цифрового фильтра. По значениям
и
можно вычислить постоянную
.
[8]
Если аналоговый
фильтр нормализован и
то
При этом
Рис. 14. Связь между частотами аналоговой и цифровой областей
Характеристики фильтра для бета-ритма
Fs1=13; % нижняя граничная частота полосы задерживания (Гц)
Fs2=32; % верхняя граничная частота полосы задерживания (Гц)
Fp1=14; % нижняя граничная частота полосы пропускания (Гц)
Fp2=30; % верхняя граничная частота полосы пропускания (Гц)
Рис.15. АЧХ И ЛАЧХ получившегося фильтра для бета-ритма, N=6
Характеристики фильтра для гамма-ритма
Частотный диапазон данного ритма довольно широкий. В полосу данного РГМ входят много сильных паразитных помех, а именно помехи на частотах в 50 и 100 Гц вызванных работой системы электроснабжения. Данные помехи происходят в силу того, что в настоящий момент напряжение в наших электросетях имеет форму синусоиды с частотой 50 Гц, которая пассивно создает электромагнитные импульсы такой же частоты. Тем самым вызывая шумы относительно высокой амплитуды.
Поэтому для выделения данного РГМ необходимо создать дополнительные режекторные фильтры (РФ), удаляющие помехи на частотах 50 и 100Гц.
Режекторный фильтр 1, для частоты 50Гц:
Fs1=49; % нижняя граничная частота полосы задерживания (Гц)
Fs2=51.5; % верхняя граничная частота полосы задерживания (Гц)
Fp1=49.5; % нижняя граничная частота полосы пропускания (Гц)
Fp2=51; % верхняя граничная частота полосы пропускания (Гц)
Переход от нормирующего ФНЧ к РФ осуществляется на основе замены в передаточной функции вида:
Рис.16. АЧХ получившегося РФ
Режекторный фильтр 2, для частоты 100Гц:
Fs1=99; % нижняя граничная частота полосы задерживания (Гц)
Fs2=101.5; % верхняя граничная частота полосы задерживания (Гц)
Fp1=99.5; % нижняя граничная частота полосы пропускания (Гц)
Fp2=101; % верхняя граничная частота полосы пропускания (Гц)
Рис.17. АЧХ получившегося РФ
Полосовой фильтр для гамма-ритма:
Fs1=31; % нижняя граничная частота полосы задерживания (Гц)
Fs2=180; % верхняя граничная частота полосы задерживания (Гц)
Fp1=33; % нижняя граничная частота полосы пропускания (Гц)
Fp2=170; % верхняя граничная частота полосы пропускания (Гц)
Рис.18. АЧХ и ЛАЧХ получившегося ПФ
Суммарный фильтр для гамма-ритма:
Рис.19 . АЧХ суммарного фильтра для гамма-ритма
Характеристики фильтра для дельта-ритма
Fs1=0.5; % нижняя граничная частота полосы задерживания (Гц)
Fs2=3.3; % верхняя граничная частота полосы задерживания (Гц)
Fp1=0.7; % нижняя граничная частота полосы пропускания (Гц)
Fp2=3; % верхняя граничная частота полосы пропускания (Гц)
Рис.20. АЧХ И ЛАЧХ получившегося фильтра для дельта-ритма N=5
Характеристики фильтра для тета-ритма
Fs1=4; % нижняя граничная частота полосы задерживания (Гц)
Fs2=7.5; % верхняя граничная частота полосы задерживания (Гц)
Fp1=4.5; % нижняя граничная частота полосы пропускания (Гц)
Fp2=7; % верхняя граничная частота полосы пропускания (Гц)
Рис.20. АЧХ И ЛАЧХ получившегося фильтра для дельта-ритма N=6
В нашей системе так же является необходимым выделение ультранизкого диапазона частот, в теории характерного для электромагнитных излучений головного мозга с частотой порядка 0.01 Гц
Для выполнения этой задачи спроектируем дополнительный фильтр низких частот (ФНЧ) с полосой пропускания до 0.5 Гц.
Дополнительный ФНЧ
Немного подкорректировав отрезок кода (2) смоделируем необходимый ФНЧ:
Fs=0.5; % нижняя граничная частота полосы задерживания (Гц)
Fp=0.01; % верхняя граничная частота полосы пропускания (Гц)
Ft=1000; % частота отсчетов (Гц)
Rs=50; % Дб
Rp=0.5; % Дб
ws=2*pi*Fs/Ft; % рад/с
wp=2*pi*Fp/Ft; % рад/с
Ws = 2*Ft*tan(ws/2);
Wp = 2*Ft*tan(wp/2);
Wst=(Ws)/(Wp);
[n, Wn] = buttord(1, Wst, Rp, Rs, 's');
[num, den] = butter(n, Wn, 's');
[BT, AT] = lp2lp(num, den, Wp);
[numd, dend] = bilinear(BT, AT, Ft);
[H, f] = freqz(numd, dend, 256, Ft);
figure(2)
subplot(211)
plot(f, abs(H)), grid
title('Magnitude response')
xlabel('Frequency, Hz ')
subplot(212)
plot(f, 20*log10(abs(H))), grid
title('Logarifmic magnitude response')
axis([0, 11000, -320, 10])
xlabel('Frequency, Hz ')
Рис.21. АЧХ И ЛАЧХ получившегося ФНЧ N=5
В результате проектирования ПФ получилось, что ритмы головного мозга будут полностью укладываться в полосах фильтров им соответствующих, так как частотный диапазон этих ритмов, полностью соответствует полосам пропускания всех спроектированных фильтров (т.е. каждый РГМ полностью входит в полосу пропускания соответствующего фильтра).
Поэтому, все спроектированные фильтры подходят для реализации системы фильтрации и разделения исходной усиленного сигнала по РГМ.