Общие основы линейной фильтрации

Фильтрация – это преобразование заданного сигнала с помощью линейного фильтра.

Рассмотрим основы линейной фильтрации на примере линейного стационарного фильтра второго порядка, передаточная функция которого имеет вид:

, (1)

где x-заданный процесс;

y – получаемый на выходе фильтра процесс;

ω₀ – частота собственных колебаний фильтра;

ξ – относительный коэффициент затухания этого фильтра.

Для контроля и графического представления передаточной функции любого динамического звена удобно использовать процедуру freqs.

В общем случае обращение к ней имеет вид:

h= freqs (b,a,w).

При этом процедура создает вектор h комплексных значений частотной характеристики W(jω) по передаточной функции W(s), заданной векторами коэффициентов ее числителя (b) и знаменателя (a), а также по заданному вектору w частоты ω. Если аргумент w не задан, то выбирается 200 отсчетов частоты, для которых вычисляется частотная характеристика.

Пример.

A=1; ξ=0.05; T₀ = 2π / ω₀=1.

Вычислим значения коэффициентов числителя и знаменателя и выведем графики АЧХ и ФЧХ:

>> T0=1;dz=0.05;

>> om0=2*pi/T0;a=1

>> A=1;

>> a1(1)=1;a1(2)=2*dz*om0;a1(3)=om0^2;b1(1)=A;

>> freqs(b1,a1)

Дискретная передаточная функция, соответствующая непрерывному фильтру (1) (которая может быть получена путем Z - преобразования разностного уравнения - дискретного аналога дифференциального уравнения), имеет вид:

, (2)

где =1+2 ξ (ω₀ Ts)+ (ω₀ Ts)²;

= -2 (1+ξ (ω₀ Ts));

=1.

В общем случае дискретная передаточная функция имеет вид:

(3)

В системе MatLab фильтрация сигнала, описываемого дискретной передаточной функцией вида осуществляется процедурой filterследующим образом:

y= filter(b, a, x)

где x – заданный вектор значений входного сигнала;

y – вектор значений выходного сигнала фильтра, получаемого вследствие фильтрации;

b – вектор коэффициентов числителя дискретной передаточной функции (1);

a – вектор коэффициентов знаменателя дискретной этой функции.

Пример

Полезный сигнал имеет синусоидальную форму с известным периодом T1 =1 c и амплитудой A1 = 0.75. сформируем этот сигнал как вектор его значений в дискретные моменты времени с дискретом Ts=0.001 c:

>> Ts=0.001;

>> t=0:Ts:20;

>> A1=0.75; T1=1;

>> Yp=A1*sin(2*pi*t/T1);

>> plot(t,Yp),grid;

>> title('P-signal');

>> xlabel('t, c');

>> ylabel('Yp')

Пусть вследствие прохождения через преобразователь и в результате измерения к полезному сигналу добавились шумы:

- более высокочастотная синусоида с периодом T2=0.2 c и амплитудой A2=10;

- белый гауссовский шум измерителя с интенсивностью Aш=5.

В результате создается такой измеренный сигнал x(t):

>> T2=0.2;A2=10; eps=pi/4;

>> Ash=5;

>> x=A1*sin(2*pi*t/T1)+A2*sin(2*pi*t/T2+eps)+Ash*randn(1,length(t));

>> plot(t,x), grid

>> title('input')

>> xlabel('t, c')

>> ylabel('x(t)')

Требуется так обработать измеренные данные x, чтобы восстановить по ним полезный процесс как можно точнее.

Т.к. частота полезного сигнала заранее известна, восстановление его можно выполнить при помощи резонансного фильтра вида (2). При этом необходимо создать такой фильтр, чтобы период его собственных колебаний Tфбыл равен периоду колебаний полезного сигнала (Tф= T1).

Чтобы после прохождения через такой фильтр амплитуда восстановленного сигнала совпадала с амплитудой полезного сигнала, нужно входной сигнал фильтра домножить на постоянную величину 2ξω₀²(поскольку при резонансе амплитуда выходного сигнала уменьшается именно во столько раз по сравнению с амплитудой входного сигнала).

Сформируем фильтр, описанный выше:

>> T1=1; Tf=T1;dz=0.05;

>> om0=2*pi/Tf;A=1;oms=om0*Ts;

>> a(1)=1+2*dz*oms+oms^2;

>> a(2)= - 2*(1+dz*oms);

>> a(3)=1;

>> b(1)=A*Ts*Ts*(2*dz*om0^2);

Пропустим сформированный процесс через него:

>> y=filter(b,a,x);

>> plot(t(10002:end),y(10002:end),t(10002:end),Yp(10002:end)),grid

>> title('procesOut(Tf=1,dz=0.05)');

>> xlabel('t, c');

>> ylabel('Y(t)')

В результате получим процесс, приведенный на рисунке.

Из рисунка можно заметить два недостатка:

- восстановленный процесс устанавливается на выходе фильтра только спустя некоторое время вследствие нулевых начальных условий самого фильтра как динамического звена;

- в установившемся режиме наблюдается значительный сдвиг фаз (π/2) между восстанавливаемым и восстановленным процессами.

Для избежания фазовых искажений полезного сигнала при его восстановлении можно воспользоваться процедурой двойной фильтрации – filtfilt, которая осуществляет обработку вектора в два приема: сначала в прямом, а затем в обратном направлении.

Результат представлен ниже.

>> y=filtfilt(b,a,x);

>> plot(t,y,t,Yp),grid;

>> title('procedure FILTFILT (Tf=1,dz=0.05)');

>> xlabel('t, c');

>> ylabel('Y(t)')