Разработка и тестирование алгоритма цифровой фильтрации

Создадим алгоритм цифровой фильтрации. Для этого воспользуемся данными, полученным с помощью метода Эйлера для периода дискретизации T_d=0,00019096 с.

b₀ = 0,1005

b₁ = -0,1005

a₁ = -1,7798

a₂ = 0,7894

Передаточная функция:

Разностное уравнение:

С помощью пакета Matlab Simulink составим структурную схему:

Рис. 31. Структурная схема алгоритма фильтрации

Этой схеме соответствует следующая функция:

Код Matlab:

function V=EuFilter(U)

b0=0.1005;

b1=-0.1005;

a1=-1.7798;

a2=0.7894;

length=size(U,2);

V=zeros(1,length);

V(1)=b0*U(1);

V(2)=b0*U(2)+b1*U(1)-a1*V(1);

for n=3:1:length

V(n)=b0*U(n)+b1*U(n-1)-a1*V(n-1)-a2*V(n-2);

end

end

В качестве тестовых сигнала используем синусоиды с разными частотами. Пропустим их через фильтр, построим графики входных и выходных сигналов, их спектров.

Чтобы легче было проанализировать полученные результаты, приведём АЧХ и ФЧХ цифровой цепи:

Код Matlab:

figure(5);

w=50:1:6580;

Wd=(b0+b1*exp(-1i*w*Td))./(1+a1*exp(-1i*w*Td)+a2*exp(-2*1i*w*Td));

subplot(2,1,1), plot(w,abs(Wd),'b');

hold on;

subplot(2,1,2), plot(w,180/pi*angle(Wd),'b');

hold on;

subplot(2,1,1), grid on, xlabel('w (Rad/s)'), title('MAGNITUDE - |H(w)|');

subplot(2,1,2), grid on, xlabel('w (Rad/s)'), title('PHASE - arg [H(w)] (deg)');

Рис. 32. АЧХ и ФЧХ фильтра

Минимальное подавление амплитуды наблюдается при частоте ω=544,5 рад/с. В этом случае выходная амплитуда будет составлять 0,478 от входной. ФЧХ пересекает линию 0º при частоте ω=577,64 рад/с. При этой частоте выходной сигнал не будет сдвинут относительно входного.

Тестирование фильтра:

1. Примем частоту ω=40 рад/с

Код Matlab:

clc;clear;

Td=0.00019096; % Период дискретизации

N=2^13; % Количество отсчётов

t=0:Td*(N-1)*Td;

w=40;

U=sin(w*t); % Тестовый сигнал

V=EuFilter(U); % Фильтрация тестового сигнала

figure(1);

plot(t(1:1024),U(1:1024),'r');

hold on;

plot(t(1:1024),V(1:1024),'g');

grid on, xlabel('t (s)'), title('Signal');

Fd=1/Td; % Частота дискретизации

SpectrU=fft(U,N);%Быстрое преобразование Фурье

SpectrU=2*SpectrU./N; % Нормировка спектра по амплитуде

SpectrU(1)=SpectrU(1)/2; % Нормировка постоянной составляющей в спектре

W=0:2*pi*Fd/N:1000;

figure(2);

subplot(2,1,1), plot(W,abs(SpectrU(1:length(W))),'r'), grid on, xlabel('w (Rad/s)'), title('MAGNITUDE');

hold on;

subplot(2,1,2), plot(W,angle(SpectrU(1:length(W)))*180/pi,'r'), grid on, xlabel('w (Rad/s)'), title('PHASE (deg)');

hold on;

SpectrV=fft(V,N); %Быстрое преобразование Фурье

SpectrV=2*SpectrV./N; % Нормировка спектра по амплитуде

SpectrV(1)=SpectrV(1)/2; % Нормировка постоянной составляющей в спектре

subplot(2,1,1), plot(W,abs(SpectrV(1:length(W))),'g');

subplot(2,1,2), plot(W,angle(SpectrV(1:length(W)))*180/pi,'g');

Рис. 33. Тестовый сигнал U=sin(40·t) (обозначен красным), и сигнал V, полученный в результате фильтрации (обозначен зелёным)

Рис. 34. Амплитудный и фазовый спектры тестового сигнала U=sin(40·t) (обозначены красным), и сигнала V, полученный в результате фильтрации (обозначены зелёным). Получены с помощью быстрого преобразования Фурье.

По спектрам определим амплитуды и фазы входного и выходного сигналов:

Амплитуда входного: A(U)≈1;

Амплитуда выходного: A(V)≈0,08;

Соотношение амплитуд:

Фаза входного: φ(U)≈-90º;

Фаза выходного: φ(V)≈-10º;

Разность фаз: φ(V)-φ(U)=-10º+90º=80º

Полученные результаты согласуются с АЧХ и ФЧХ фильтра для этой частоты. Она находится левее основной полосы пропускания.

2. Примем частоту ω=545 рад/с

Рис. 35. Тестовый сигнал U=sin(545·t) (обозначен красным), и сигнал V, полученный в результате фильтрации (обозначен зелёным)

Рис. 36. Амплитудный и фазовый спектры тестового сигнала U=sin(545·t) (обозначены красным), и сигнала V, полученный в результате фильтрации (обозначены зелёным).

По спектрам определим амплитуды и фазы входного и выходного сигналов:

Амплитуда входного: A(U)≈1;

Амплитуда выходного: A(V)≈0,473;

Соотношение амплитуд:

Фаза входного: φ(U)≈-90º;

Фаза выходного: φ(V)≈-87º;

Разность фаз: φ(V)-φ(U)=-87º+90º=3º

Полученные результаты согласуются с АЧХ и ФЧХ фильтра для этой частоты. Она находится в середине основной полосы пропускания.

3. Примем частоту ω=1000 рад/с

Рис. 37. Тестовый сигнал U=sin(1000·t) (обозначен красным), и сигнал V, полученный в результате фильтрации (обозначен зелёным)

Рис. 38. Амплитудный и фазовый спектры тестового сигнала U=sin(1000·t) (обозначены красным), и сигнала V, полученный в результате фильтрации (обозначены зелёным).

По спектрам определим амплитуды и фазы входного и выходного сигналов:

Амплитуда входного: A(U)≈1;

Амплитуда выходного: A(V)≈0,4;

Соотношение амплитуд:

Фаза входного: φ(U)≈-90º;

Фаза выходного: φ(V)≈-114º;

Разность фаз: φ(V)-φ(U)=-114º+90º=-24º

Полученные результаты согласуются с АЧХ и ФЧХ фильтра для этой частоты. Она находится в правой части основной полосы пропускания.

Подадим на вход фильтра равномерный белый шум:

Код Matlab:

U=rand(1,N);