96958

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

Курсовая

Коммуникация, связь, радиоэлектроника и цифровые приборы

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

Русский

2015-10-12

1.98 MB

2 чел.

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

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

высшего профессионального образования

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

электроники и автоматики"

МГТУМИРЭА

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

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

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

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

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

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

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

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

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

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

профессор, д.т.н., профессор

Асанов А.З.

Работа представлена к защите

«__»_______2015 г.

(подпись студента)

«Допущен к защите»

«__»_______2015 г.

(подпись руководителя)

Москва 2015

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

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

высшего профессионального образования

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

электроники и автоматики"

МГТУМИРЭА

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

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

Утверждаю

Заведующий

кафедрой_____________

«____» __________201___ г.

ЗАДАНИЕ

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

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

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

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

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

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

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

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

работу выдал

«___»______201__г.

Подпись руководителя

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

работы

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

работу получил

«___»______201__г.

Подпись студента –

исполнителя

Ф.И.О. студента -

исполнителя

работы

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

№ этапа

Этап курсовой работы выполнил и представил результаты руководителю проекта (работы),

дата и подпись исполнителя

Работу по этапу курсовой  работы принял на рассмотрение,

дата и подпись руководителя

Рекомендации и замечания по
этапу курсовой

работы
выдал
исполнителю,

дата и подпись руководителя

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

1

2

3

4

5

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

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

высшего профессионального образования

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

электроники и автоматики"

МГТУМИРЭА

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

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

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

от ________________201__г.  №__________

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

______________________________________________________________________

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

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

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

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

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

____________________________________________________________________________________________

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

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

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

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

1 _____________________________________________________________________________

2. _____________________________________________________________________________

3. _____________________________________________________________________________

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

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

____________________________________________

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


Оглавление

[1] Федеральное государственное бюджетное образовательное учреждение

[2] высшего профессионального образования

[3] "Московский государственный технический университет радиотехники,

[4] электроники и автоматики"

[5] МГТУМИРЭА

[6] Федеральное государственное бюджетное образовательное учреждение

[7] высшего профессионального образования

[8] "Московский государственный технический университет радиотехники,

[9] электроники и автоматики"

[10] МГТУМИРЭА

[11] Федеральное государственное бюджетное образовательное учреждение

[12] высшего профессионального образования

[13] "Московский государственный технический университет радиотехники,

[14] электроники и автоматики"

[15] МГТУМИРЭА

[16] Оглавление

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

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

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

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

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


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

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

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

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

R1=1000 Ом

R2=4000 Ом

С1=1,5 мкФ

C2=0,5 мкФ

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

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

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

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

Далее найдём потенциал в точке 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

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

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

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

C0=1,5

C1=-0,5

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

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

Код в Matlab:

figure(2);

impulse(Wanalog,'g');

grid on;

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

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

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

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

ωmax=6580,668 рад/с

fmax= =1047,346 Гц

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

fd=5fmax=5236,730 Гц

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

Td=  =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];

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

b0 =  0,1005

b1 = -0,1005

a1 =  -1,7798

a2 = 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. Карта нулей и полюсов цифрового фильтра, полученного методом Эйлера при Td=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.  ЛАЧХ аналогового фильтра (обозначена зелёным) и соответствующего ему цифрового фильтра, полученного методом Эйлера при Td=0,00019096 с (обозначена синим)

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

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

Как видно из полученных графиков (рисунки 7 и 9), частотные характеристики цифрового фильтра по форме почти точно повторяют характеристики аналогового прототипа на участке [0; ], то есть [0; 1,645·104] рад/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.  Импульсная переходная функция цифрового фильтра, полученного методом Эйлера при Td=0,00019096 с

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

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

fd=10fmax=10473,460 Гц

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

Td=  =0,00009548с

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

b0 =  0,0563

b1 = -0,0563

a1 =  -1,8820

a2 = 0,8847

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

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

Нули:

z=0; z=1

Полюса:

z=0,9692; z= 0,9128

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

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

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

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

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

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

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

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

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

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

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

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

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

Td=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];

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

b0 =  0,0563

b1 = 0

b2 = -0,0563

a1 =  -1,7640

a2 = 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. Карта нулей и полюсов цифрового фильтра, полученного методом билинейного преобразования при Td=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.  ЛАЧХ аналогового фильтра (обозначена зелёным) и соответствующего ему цифрового фильтра, полученного методом билинейного преобразования при Td=0,00019096 с (обозначена синим)

Рис. 17.  ЛАЧХ и ЛФЧХ цифрового фильтра, полученного методом билинейного преобразования при Td=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.  Импульсная переходная функция цифрового фильтра, полученного методом билинейного преобразования при Td=0,00019096 с

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

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

Td=0,00009548с

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

b0 =  0,0299

b1 = 0

b2 = -0,0299

a1 =  -1,8775

a2 = 0,8804

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

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

Нули:

z=-1; z=1

Полюса:

z=0,9687; z= 0,9089

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Td=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];

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

b0 =  0,1273

b1 = -0,1266

a1 =  -1,7645

a2 = 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. Карта нулей и полюсов цифрового фильтра, полученного методом инвариантной импульсной характеристики при Td=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.  ЛАЧХ аналогового фильтра (обозначена зелёным) и соответствующего ему цифрового фильтра, полученного методом инвариантной импульсной характеристики при Td=0,00019096 с (обозначена синим)

Рис. 25.  ЛАЧХ и ЛФЧХ цифрового фильтра, полученного методом инвариантной импульсной характеристики при Td=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.  Импульсная переходная функция цифрового фильтра, полученного методом инвариантной импульсной характеристики при Td=0,00019096 с

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

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

Td=0,00009548с

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

b0 =  0,0637

b1 = -0,0636

a1 =  -1,8776

a2 = 0,8805

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

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

Нули:

z=0; z=0,9985

Полюса:

z=0,9687; z= 0,9089

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

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

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

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

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

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

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

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

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

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

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


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

  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º

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

  1.  Примем частоту ω=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º

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

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

PAGE   \* MERGEFORMAT43


 

А также другие работы, которые могут Вас заинтересовать

68976. Умовний оператор. Оператор вибору. Цикли 38 KB
  Виконання тіла оператора-перемикача switch починається з вибраного таким чином оператора і продовжується до кінця тіла або до тих пір, поки який-небудь оператор не передасть управління за межі тіла. Оператор, наступний за ключовим словом default, виконується, якщо жодна з...
68977. Одновимірні та багатовимірні масиви 30 KB
  Відповідно до синтаксису Сі в мові існують тільки одновимірні масиви, проте елементами одновимірного масиву, у свою чергу, можуть бути масиви. Тому двовимірний масив визначається як масив масивів. Таким чином, в прикладі визначений масив Z з 13 елементів-масивів, кожний з яких...
68978. Вказівники. Функції динамічного розподілу пам’яті 37 KB
  Кожна змінна в програмі - це об’єкт, який має ім’я і значення. За ім’ям можна звернутися до змінної і отримати (а потім, наприклад, надрукувати) її значення. Щоб отримати адресу в явному вигляді, в мові Сі застосовують унарну операцію. Вираз Е дозволяє отримати адресу ділянки пам’яті, виділеної на машинному рівні для змінної Е.
68979. Функції, їх параметри. Рекурсія. Прототипи функцій 35.5 KB
  Визначення функції Опис функції та її тип Рекурсивні функції Визначення функції. Синонімами цього іншого поняття в інших мовах програмування є процедури підпрограми підпрограми-функції процедури-функції. Всі функції в мові Сі мають рекомендуємий стандартами мови єдиний формат...
68980. Структури, об’єднання 36.5 KB
  Структура - це з’єднане в єдине ціле безліч поіменованих елементів (компонентів) даних. На відміну від масиву, який завжди складається з однотипних елементів, компоненти структури можуть бути різних типів і всі повинні мати різні імена.
68981. Рекурсивні функції і процедури, параметри-процедури 30 KB
  Тобто це є визначенням функції через цю саму функцію, У мові Паскаль рекурсивний опис функції полягає в тому, що в тілі такої функції міститься звертання до цієї ж функції. Наведемо рекурсивний опис функції п...
68982. Файли, робота з файлами 41 KB
  План заняття: Організація файлів Робота з файлами Підготовчі та завершальні операції Операції уведеннявиведення Пересування по файлу Організація файлів Є багато задач коли кількість компонентів певного типу будьякого з відомих уже нам наперед визначити неможливо то її визначають у процесі виконання програми.
68983. Текстові файли 36.5 KB
  В кінці кожного рядка є символ кінець рядка внутрішнє відображення якого залежить від реалізації. Звичайно кінець рядка це комбінація коду переведення каретки символ 13 за яким може бути код переведення рядка символ 10. Для програмування переважно немає потреби знати коди символів...
68984. Модулі. Модуль і його структура 49.5 KB
  Модуль - це сукупність сталих, типів даних, змінних, процедур і функцій, які можна використати у програмі або в іншому модулі. Сам модуль не є виконуваною програмою. Модульний підхід до проектування дає змогу розділити програму на частини, які компілюють окремо.