- •5. Мониторинг процесса выполнения курсовой работы
- •Протокол заседания комиссии по защите курсовой работы
- •Исследование характеристик аналоговой цепи
- •Расчёт цифровой цепи методом Эйлера.
- •Расчёт цифровой цепи методом билинейного преобразования.
- •Расчёт цифровой цепи методом инвариантной импульсной характеристики.
- •Разработка и тестирование алгоритма цифровой фильтрации.
Разработка и тестирование алгоритма цифровой фильтрации.
Создадим алгоритм цифровой фильтрации. Для этого воспользуемся данными, полученным с помощью метода Эйлера для периода дискретизации Td=0,00019096 с.
b0 = 0,1005
b1 = -0,1005
a1 = -1,7798
a2 = 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 рад/с. При этой частоте выходной сигнал не будет сдвинут относительно входного.
Тестирование фильтра:
Примем частоту ω=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º
Полученные результаты согласуются с АЧХ и ФЧХ фильтра для этой частоты. Она находится левее основной полосы пропускания.
Примем частоту ω=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º
Полученные результаты согласуются с АЧХ и ФЧХ фильтра для этой частоты. Она находится в середине основной полосы пропускания.
Примем частоту ω=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);
Рис. 39. Белый шум до фильтрации (обозначен красным) и после (обозначен зелёным)
Рис. 40. Спектр белого шума до фильтрации (обозначен красным) и после (обозначен зелёным)
Подадим на вход сигнал, состоящий из трёх гармоник на частотах 40 рад/с, 545 рад/с 5000 рад/с и нормального белого шума.
Код Matlab:
clc;clear;
Td=0.00019096; %Период дискретизации
N=2^13; %Количество отсчётов
t=0:Td:(N-1)*Td;
w1=40;
w2=545;
w3=5000;
U=sin(w1*t)+sin(w2*t)+sin(w3*t)+randn(1,N);
V=EuFilter(U);
figure(1);
plot(t(1:1024),U(1:1024),'r');
hold on;
plot(t(1:1024),V(1:1024),'g','LineWidth',2);
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:5500;
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','LineWidth',2);
subplot(2,1,2), plot(W,angle(SpectrV(1:length(W)))*180/pi,'g','LineWidth',2);
Рис. 41. Тестовый сигнал U(t)= sin(40·t) + sin(545·t) +sin(5000·t)+n(t) (обозначен красным), и сигнал V, полученный в результате фильтрации (обозначен зелёным)
Рис. 42. Амплитудный спектр тестового сигнала U(t) до фильтрации (обозначен красным) и после (обозначен зелёным)
После прохождения фильтра существенное ослабление претерпела первая гармоника (с частотой 40 рад/с), находящаяся левее основной полосы пропускания фильтра. Также существенно подавлена третья гармоника (с частотой 5000 рад/с), находящаяся правее основной полосы. Кроме того, на выходе наблюдается меньшее отношение шума к полезному сигналу.