Построение и исследование характеристик цифрового эквивалента аналогового фильтра

Размещено на http://www.allbest.ru/

МИНОБРНАУКИ РОССИИ

Федеральное государственное бюджетное образовательное учреждение высшего профессионального образования

"Московский государственный технический университет радиотехники, электроники и автоматики"

МГТУМИРЭА

факультет кибернетики

Кафедра автоматических систем

КУРСОВАЯ РАБОТА

по дисциплине

«Цифровая обработка сигналов»

Тема курсовой работы

«Построение и исследование характеристик цифрового эквивалента аналогового фильтра»

Студент группы КУБ-1-12

Барашков А.А. (шифр 121290)

Руководитель курсовой работы

профессор, д.т.н., профессор Асанов А.З.

Москва 2015

Задание

на выполнение курсовой работы

по дисциплине «Цифровая обработка сигналов »

Студент Барашков Алексей Андреевич Группа КУБ-1-12

1. Тема « Построение и исследование характеристик цифрового эквивалента аналогового фильтра».

2. Исходные данные:

R1=1000 Ом, R2=4000 Ом, С1=1,5 мкФ, C2=0,5 мкФ.

3. Перечень вопросов, подлежащих разработке, и обязательного графического материала:

4. Срок представления к защите курсовой работы:до «___» _______2015 г.

Задание на курсовую

работу выдал «___»______201__г. Подпись руководителя

Ф.И.О. руководителя

работы

Задание на курсовую

работу получил «___»______201__г. Подпись студента -

исполнителя Ф.И.О. студента -

исполнителя

работы

5. Мониторинг процесса выполнения курсовой работы

№ этапа	Этап курсовой работы выполнил и представил результаты руководителю проекта (работы), дата и подпись исполнителя	Работу по этапу курсовой работы принял на рассмотрение, дата и подпись руководителя	Рекомендации и замечания по этапу курсовой работы выдал исполнителю, дата и подпись руководителя	Комментарии руководителя курсовой проекта работы
1
2
3
4
5

Протокол заседания комиссии по защите курсовой работы

от ________________201__г. №__________

Состав комиссии:___________________________________________

_____________________________________________________________

(должность, ученая степень, ученое звание)

Утверждена распоряжением заведующего кафедрой ________________

(наименование кафедры)

от «___» ________ 201_ г. №________.

Слушали защиту курсовой работы __________________________

_____________________________________________________________

по дисциплине ____________________________________________

студента группы______________________________________________

(группа) (Ф.И.О.)

Во время защиты курсовой работы были заданы следующие вопросы:

1 ___________________________________________________________

2. __________________________________________________________

3. ___________________________________________________________

Итоговая (комплексная) оценка выполнения и защиты курсовой работы_____________________

Члены комиссии ____________________________________________

____________________________________________

(подпись) (Ф.И.О.)

Содержание

Исследование характеристик аналоговой цепи

Расчёт цифровой цепи методом Эйлера

Расчёт цифровой цепи методом билинейного преобразования

Расчёт цифровой цепи методом инвариантной импульсной характеристики

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

Исследование характеристик аналоговой цепи

Дана аналоговая пассивная цепь:

Рис. 1. Схема аналоговой цепи

Для цепи заданы параметры:

R₁=1000 Ом

R₂=4000 Ом

С₁=1,5 мкФ

C₂=0,5 мкФ

Воспользуемся символическим методом расчёта электрических цепей.

Рис. 2. Схема для расчёта аналоговой цепи

Применим метод «двух узлов». Узел a - заземлим, его потенциал будет равен нулю. Найдём потенциал узла c.

Ток, протекающий через резистор R₂:

Далее найдём потенциал в точке b. Так как , а ,

Получим частотную передаточную функцию цепи:

Код в Matlab:

clc;clear;

R1=1000;

R2=4000;

C1=1.5*10^(-6);

C2=0.5*10^(-6);

Wa0=(R2*C2)^2;

Wa1=R2*C2;

Wb0=C1*R1*(R2*C2)^2;

Wb1=C2*R2*(C2*R2+2*C1*R1+C2*R1);

Wb2=2*C2*R2+C1*R1+C2*R1;

Wanalog=tf([Wa0 Wa1 0], [Wb0 Wb1 Wb2 1])

Подставив значения параметров R и C, получим:

Так как преобразования Фурье и Лапласа имеют схожий вид, можно выполнить замену p=jщ. Тогда передаточная функция аналоговой цепи:

Найдём для этой передаточной функции нули и полюса.

Код в Matlab:

disp('Полюса:');

disp(pole(Wanalog));

disp('Нули:');

disp(zero(Wanalog));

Нули:

p=0; p=-500

Полюса:

p=-1000; p=-500; p=-1000/3

Исходя из найденных значений, можем записать передаточную функцию в виде:

Этой передаточной функции соответствует дифференциальное уравнение:

Получим логарифмические амплитудно-частотные и фазо-частотные характеристики для этого фильтра:

Амплитудно-частотная характеристика:

Фазо-частотная характеристика:

Построим графики этих характеристик:

Код в Matlab:

figure(1);

bode(Wanalog,'g');

grid on;

Рис. 3. Логарифмические частотные характеристики аналогового фильтра

Для проверки предыдущих вычислений можно сравнить полученную выше ЛАЧХ с ЛАЧХ, полученной при непосредственном моделировании электрической цепи в пакете Electronics Workbench:

Рис. 4. Моделирование электрической цепи в пакете Electronics Workbench

Получим импульсную переходную функцию этого фильтра.

Для этого разложим передаточную функцию на простые дроби:

Найдём коэффициенты:

C₀=1,5

C₁=-0,5

С помощью обратного преобразования Лапласа, найдём ИПФ:

Построим в системе Matlab график ИПФ:

Код в Matlab:

figure(2);

impulse(Wanalog,'g');

grid on;

Рис. 5. График импульсной переходной функции аналогового фильтра.

Расчёт цифровой цепи методом Эйлера.

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

По ЛАЧХ аналогового фильтра найдём частоту, при которой происходит ослабление амплитуды в 10 раз (-20дБ). Примем её за ширину спектра.

щ_max=6580,668 рад/с

f_max= =1047,346 Гц

Возьмём частоту дискретизации в 5 раз больше ширины спектра:

f_d=5f_max=5236,730 Гц

Тогда период дискретизации:

T_d= =0,00019096 с

Выбрав период дискретизации, можно найти коэффициенты передаточной функции и разностного уравнения:

Код Matlab:

clc;clear;

Td=0.00019096;

%Период дискретизации

A=2/3*10^3;

%Полюса аналоговой передаточной функции

p0=1000;

p1=1000/3;

u=(1+Td*(p0+p1+Td*p0*p1));

%Коэффициенты передаточной функции цифровой цепи

b0=A*Td/u

b1=-A*Td/u

a1=(-2-Td*(p0+p1))/u

a2=1/u

b=[b0 b1];

a=[1 a1 a2];

Полученные коэффициенты:

b₀ = 0,1005

b₁ = -0,1005

a₁ = -1,7798

a₂ = 0,7894

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

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

Найдём нули и полюса этой передаточной функции:

Код Matlab:

[q,p]=tf2zpk(b,a);

disp('Нули');

disp(q);

disp('Полюса');

disp(p);

figure(4);

zplane(b,a);

Нули:

z=0; z=1

Полюса:

z=0,9402; z= 0,8397

Рис. 6. Карта нулей и полюсов цифрового фильтра, полученного методом Эйлера при T_d=0,00019096 с

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

Амплитудно-частотная характеристика:

Фазо-частотная характеристика:

Код Matlab:

w=logspace(1,5,10000);

figure(1);

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

% Построение графиков ЛАЧХ и ЛФЧХ цифровой цепи

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

hold on;

subplot(2,1,2), semilogx(w,180/pi*angle(Wd),'b'), grid on, xlabel('w (Rad/s)'), title('PHASE - arg [H(w)] (deg)');

hold on;

figure(2);

% Построение графиков ЛАЧХ аналоговой и цифровой цепи на одном полотне Wa=A*1i*w./(1i*w+p0)./(1i*w+p1);

loglog(w,abs(Wd),'b');

hold on;

loglog(w,abs(Wa),'g');

hold on;

grid on; xlabel('w (Rad/s)'); title('MAGNITUDE - |H(w)|'); axis([0 10^5 10^(-2) 1]);

Рис. 7. ЛАЧХ аналогового фильтра (обозначена зелёным) и соответствующего ему цифрового фильтра, полученного методом Эйлера при T_d=0,00019096 с (обозначена синим)

Рис. 8. ЛАЧХ и ЛФЧХ цифрового фильтра, полученного методом Эйлера при T_d=0,00019096 с

Рис. 9. Амплитудно-частотная характеристика (не в логарифмическом масштабе) аналогового фильтра (обозначена зелёным) и соответствующего ему цифрового фильтра, полученного методом Эйлера при T_d=0,00019096 с (обозначена синим)

Как видно из полученных графиков (рисунки 7 и 9), частотные характеристики цифрового фильтра по форме почти точно повторяют характеристики аналогового прототипа на участке [0; ], то есть [0; 1,645·10⁴] рад/c. Далее идёт периодическое повторение характеристики цифрового фильтра. аналоговый цифровой цепь фильтрация

Импульсная переходная функция цифрового фильтра может быть получена за счёт обратного Z-преобразования его передаточной функции. Мы получим её с помощью пакета Matlab.

Код Matlab:

N=0.012/Td;

n=0:(N-1);

h=impz(b,a,N);

figure(3);

title('Impulse Response h(n*Td) - impz');

hold on;

xlabel('n');

ylabel('h(n*Td)');

plot(n,h,'b');

grid on;

Рис. 10. Импульсная переходная функция цифрового фильтра, полученного методом Эйлера при T_d=0,00019096 с

По рисунку 10 видно, что ИПФ цифрового фильтра очень близка по форме к ИПФ аналогового и отличается, в первую очередь, на величину вещественного коэффициента.

Теперь возьмём частоту дискретизации в 10 раз больше ширины спектра:

f_d=10f_max=10473,460 Гц

Тогда период дискретизации:

T_d= =0,00009548с

Найдём коэффициенты передаточной функции и разностного уравнения:

b₀ = 0,0563

b₁ = -0,0563

a₁ = -1,8820

a₂ = 0,8847

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

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

Нули:

z=0; z=1

Полюса:

z=0,9692; z= 0,9128

Рис. 11. Карта нулей и полюсов цифрового фильтра, полученного методом Эйлера при T_d=0,00009548 с

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

Амплитудно-частотная характеристика:

Фазо-частотная характеристика:

Рис. 12. ЛАЧХ аналогового фильтра (обозначена зелёным) и соответствующего ему цифрового фильтра, полученного методом Эйлера при T_d=0,00009548 с (обозначена синим)

Рис. 13. ЛАЧХ и ЛФЧХ цифрового фильтра, полученного методом Эйлера при T_d=0,00009548 с

В результате сравнения графиков на рисунках 7 и 12, 8 и 13 можно сделать вывод, что при уменьшении периода дискретизации различия между частотными характеристиками цифрового фильтра и его аналогового прообраза становятся несколько меньше. Также у частотных характеристик цифрового фильтра увеличивается период.

Импульсная переходная функция:

Рис. 14. Импульсная переходная функция цифрового фильтра, полученного методом Эйлера при T_d=0,00009548 с

Из сравнения ИПФ при разных периодах дискретизации (рис. 10 и рис. 14) видно, что при уменьшении T_d в несколько раз, во столько же раз уменьшается амплитуда ИПФ.

Расчёт цифровой цепи методом билинейного преобразования.

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

Примем период дискретизации:

T_d=0,00019096 с

Выбрав период дискретизации, можно найти коэффициенты передаточной функции и разностного уравнения:

Код Matlab:

clc;clear;

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

A=2/3*10^3;

%Полюса аналоговой передаточной функции

p0=1000;

p1=1000/3;

u=(2/Td+p0)*(2/Td+p1);

%Коэффициенты передаточной функции цифровой цепи

b0=2*A/Td/u

b1=0

b2=-2*A/Td/u

a1=(-8/Td^2+2*p0*p1)/u

a2=(-2/Td+p0)*(-2/Td+p1)/u

b=[b0 b1 b2];

a=[1 a1 a2];

Полученные коэффициенты:

b₀ = 0,0563

b₁ = 0

b₂ = -0,0563

a₁ = -1,7640

a₂ = 0,7747

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

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

Найдём нули и полюса этой передаточной функции:

Код Matlab:

[q,p]=tf2zpk(b,a);

disp('Нули');

disp(q);

disp('Полюса');

disp(p);

figure(4);

zplane(b,a);

Нули:

z=-1; z=1

Полюса:

z=0,9383; z= 0,8257

Рис. 15. Карта нулей и полюсов цифрового фильтра, полученного методом билинейного преобразования при T_d=0,00019096 с

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

Амплитудно-частотная характеристика:

Фазо-частотная характеристика:

Код Matlab:

w=logspace(1,5,10000);

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

% Построение графиков ЛАЧХ и ЛФЧХ цифровой цепи

figure(1);

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

hold on;

subplot(2,1,2), semilogx(w,180/pi*angle(Wd),'b'), grid on, xlabel('w (Rad/s)'), title('PHASE - arg [H(w)] (deg)');

hold on;

% Построение графиков ЛАЧХ аналоговой и цифровой цепи на одном полотне

figure(2);

Wa=A*1i*w./(1i*w+p0)./(1i*w+p1);

loglog(w,abs(Wd),'b');

hold on;

loglog(w,abs(Wa),'g');

hold on;

grid on; xlabel('w (Rad/s)'); title('MAGNITUDE - |H(w)|'); axis([0 10^5 10^(-2) 1]);

Рис. 16. ЛАЧХ аналогового фильтра (обозначена зелёным) и соответствующего ему цифрового фильтра, полученного методом билинейного преобразования при T_d=0,00019096 с (обозначена синим)

Рис. 17. ЛАЧХ и ЛФЧХ цифрового фильтра, полученного методом билинейного преобразования при T_d=0,00019096 с

Полученная АЧХ цифрового фильтра, синтезированного методом билинейного преобразования, больше соответствует АЧХ аналогового прототипа в области частот от 10 до 10^4 рад/c, чем методом Эйлера при том же периоде дискретизации.

Импульсная переходная функция цифрового фильтра может быть получена за счёт обратного Z-преобразования его передаточной функции. Мы получим её с помощью пакета Matlab.

Код Matlab:

N=0.012/Td;

n=0:(N-1);

h=impz(b,a,N);

figure(3);

title('Impulse Response h(n*Td) - impz');

hold on;

xlabel('n');

ylabel('h(n*Td)');

plot(n,h,'b');

grid on;

Рис. 18. Импульсная переходная функция цифрового фильтра, полученного методом билинейного преобразования при T_d=0,00019096 с

ИПФ фильтра, полученного методом БП, несколько отличается по форме от ИПФ прототипа для начальных отсчётов.

Теперь возьмём другой период дискретизации

T_d=0,00009548с

Найдём коэффициенты передаточной функции и разностного уравнения:

b₀ = 0,0299

b₁ = 0

b₂ = -0,0299

a₁ = -1,8775

a₂ = 0,8804

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

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

Нули:

z=-1; z=1

Полюса:

z=0,9687; z= 0,9089

Рис. 19. Карта нулей и полюсов цифрового фильтра, полученного методом билинейного преобразования при T_d=0,00009548 с

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

Амплитудно-частотная характеристика:

Фазо-частотная характеристика:



Рис. 20. ЛАЧХ аналогового фильтра (обозначена зелёным) и соответствующего ему цифрового фильтра, полученного методом билинейного преобразования при T_d=0,00009548 с (обозначена синим)

Рис. 21. ЛАЧХ и ЛФЧХ цифрового фильтра, полученного методом билинейного преобразования при T_d=0,00009548 с

Импульсная переходная функция:

Рис. 22. Импульсная переходная функция цифрового фильтра, полученного методом билинейного преобразования при T_d=0,00009548 с

При уменьшении периода дискретизации у фильтра, полученного методом билинейного преобразования, характеристики изменились так же, как и у фильтра, полученного методом Эйлера.

Расчёт цифровой цепи методом инвариантной импульсной характеристики

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

В этом методе синтеза ИПФ ЦФ определяется через дискретные значения ИПФ аналогового прототипа:

Восстановим передаточную функцию по ИПФ:

Примем период дискретизации:

T_d=0,00019096 с

Выбрав период дискретизации, можно найти коэффициенты передаточной функции и разностного уравнения:

Код Matlab:

clc;clear;

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

A=2/3*10^3;

C0=1.5;

C1=-0.5;

%Полюса аналоговой передаточной функции

p0=1000;

p1=1000/3;

%Коэффициенты передаточной функции цифровой цепи

b0=A*Td*(C0+C1)

b1=A*Td*(-C0*exp(-p1*Td)-C1*exp(-p0*Td))

a1=-exp(-p1*Td)-exp(-p0*Td)

a2=exp(-(p0+p1)*Td)

b=[b0 b1];

a=[1 a1 a2];

Полученные коэффициенты:

b₀ = 0,1273

b₁ = -0,1266

a₁ = -1,7645

a₂ = 0,7752

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

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

Найдём нули и полюса этой передаточной функции:

Код Matlab:

[q,p]=tf2zpk(b,a);

disp('Нули');

disp(q);

disp('Полюса');

disp(p);

figure(4);

zplane(b,a);

Нули:

z=0; z=0,9944

Полюса:

z=0,9383; z= 0,8262

Рис. 23. Карта нулей и полюсов цифрового фильтра, полученного методом инвариантной импульсной характеристики при T_d=0,00019096 с

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

Амплитудно-частотная характеристика:

Фазо-частотная характеристика:

Код Matlab:

w=logspace(1,5,10000);;

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

% Построение графиков ЛАЧХ и ЛФЧХ цифровой цепи figure(1)

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

hold on;

subplot(2,1,2), semilogx(w,180/pi*angle(Wd),'b'), grid on, xlabel('w (Rad/s)'), title('PHASE - arg [H(w)] (deg)');

hold on;

% Построение графиков ЛАЧХ аналоговой и цифровой цепи на одном полотне

figure(2);

Wa=A*1i*w./(1i*w+p0)./(1i*w+p1);

loglog(w,abs(Wd),'b');

hold on;

loglog(w,abs(Wa),'g');

hold on;

grid on; xlabel('w (Rad/s)'); title('MAGNITUDE - |H(w)|'); axis([0 10^5 10^(-2) 1]);

Рис. 24. ЛАЧХ аналогового фильтра (обозначена зелёным) и соответствующего ему цифрового фильтра, полученного методом инвариантной импульсной характеристики при T_d=0,00019096 с (обозначена синим)

Рис. 25. ЛАЧХ и ЛФЧХ цифрового фильтра, полученного методом инвариантной импульсной характеристики при T_d=0,00019096 с

АЧХ цифрового фильтра, синтезированного методом инвариантной импульсной характеристики, имеет существенно большие различия с АЧХ аналогового прототипа, чем у фильтров, полученных методами Эйлера и билинейного преобразования.

Импульсная переходная функция цифрового фильтра может быть получена за счёт обратного Z-преобразования его передаточной функции. Мы получим её с помощью пакета Matlab.

Код Matlab:

N=0.012/Td;

n=0:(N-1);

h=impz(b,a,N);

figure(3);

title('Impulse Response h(n*Td) - impz');

hold on;

xlabel('n');

ylabel('h(n*Td)');

plot(n,h,'b');

grid on;

Рис. 26. Импульсная переходная функция цифрового фильтра, полученного методом инвариантной импульсной характеристики при T_d=0,00019096 с

ИПФ цифрового фильтра, полученного методом ИИХ отличается от ИПФ аналогового прототипа на вещественный множитель, равный периоду дискретизации.

Теперь возьмём другой период дискретизации

T_d=0,00009548с

Найдём коэффициенты передаточной функции и разностного уравнения:

b₀ = 0,0637

b₁ = -0,0636

a₁ = -1,8776

a₂ = 0,8805

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

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

Нули:

z=0; z=0,9985

Полюса:

z=0,9687; z= 0,9089

Рис. 27. Карта нулей и полюсов цифрового фильтра, полученного методом инвариантной импульсной характеристики при T_d=0,00009548 с

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

.

Амплитудно-частотная характеристика:

Фазо-частотная характеристика:

Рис. 28. ЛАЧХ аналогового фильтра (обозначена зелёным) и соответствующего ему цифрового фильтра, полученного методом инвариантной импульсной характеристики при T_d=0,00009548 с (обозначена синим)

Рис. 29. ЛАЧХ и ЛФЧХ цифрового фильтра, полученного методом инвариантной импульсной характеристики при T_d=0,00009548 с

Импульсная переходная функция:

Рис. 30. Импульсная переходная функция цифрового фильтра, полученного методом инвариантной импульсной характеристики при T_d=0,00009548 с

При уменьшении периода дискретизации у фильтра, полученного методом билинейного преобразования, характеристики изменились так же, как и у фильтра, полученного методом Эйлера.

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

Создадим алгоритм цифровой фильтрации. Для этого воспользуемся данными, полученным с помощью метода Эйлера для периода дискретизации 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);

Рис. 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 рад/с), находящаяся правее основной полосы. Кроме того, на выходе наблюдается меньшее отношение шума к полезному сигналу.

Размещено на Allbest.ru