close

Вход

Забыли?

вход по аккаунту

?

Курсовая работа

код для вставкиСкачать
ФГБОУ ВПО «Чувашский государственный университет
имени И. Н. Ульянова»
Кафедра промышленной электроники
Курсовая работа
по теории сигналов
«Синтез цифровых фильтров»
Вариант 15
Выполнил:
ст. гр. РТЭ 31-09
Шурбин С. Н.
Проверил:
Малинин Г.В.
Чебоксары 2012
Содержание
1. Аннотация ……………………………………………………………....…3
2. Задание на курсовую работу ……………………………………………..4
3. Проектирование КИХ-фильтра …………………………………………..6
4. Проектирование БИХ-фильтра ………………………………..………...20
5. Проектирование и анализ фильтров в графической среде FDATool.....33
6. Листинг MatLab-программы …………………………………………….39
7. Заключение …………………………………………………………...…..47
8. Список использованной литературы ………………………....................48
2
Аннотация
В данной курсовой работе производится расчет цифровых полосовых
фильтров с конечной и бесконечной импульсными характеристиками. Для
проверки расчета используются специализированные MatLab-функции, а так
же графическая среда для расчета и анализа фильтров FDATool.
3
Задание на курсовую работу
Таблица 1
АФ-прототип
Тип
фильтр АЧХ
а
ПФ
Б
ЦФ
Rp ,
Rs ,
дБ
дБ
1,5
40
Fp , кГц
Fs , кГц
5,0…9,0 3,0…11,0
Частота
Метод дискретизаци
синтеза
и f д , кГц
ИП
50
Используемые сокращения:
Б – фильтр Баттерворта
R p – допустимый уровень пульсаций АЧХ в полосе пропускания
Rs – минимально необходимое затухание в полосе задержки
Fp , Fs – границы полос пропускания и задержки
ИП – инвариантная импульсная характеристика
1. Проектирование БИХ – фильтра
В соответствии с таблицей 1 необходимо:
1) С помощью системы MatLab синтезировать АФ-прототип. Записать
передаточную функцию АФ. Рассчитать и построить частотные
характеристики АФ-прототипа. Сравнить АЧХ АФ с АЧХ идеального
фильтра.
2) Синтезировать рекурсивный ЦФ по аналоговому прототипу. Синтез
ЦФ провести без использования специализированных функций MatLab
синтеза ЦФ и с их использованием. Этапы синтеза отразить в работе.
3) Записать передаточную функцию ЦФ. Построить диаграмму полюсов
и нулей.
4) Составить разностное уравнение ЦФ. Решить разностное уравнение
методом итерационной процедуры для управляющего сигнала типа
единичный скачок. Построить кривую переходного процесса. С помощью
системы MatLab построить импульсную и переходную характеристики
фильтра и сигнала на выходе для произвольного входа.
5) Определит устойчивость фильтра.
6) Составить каноническую структурную схему фильтра. Разложением
передаточной функции ЦФ на произведение/сумму звеньев второго (первого)
порядка реализовать структурную схему фильтра в последовательной или
параллельной форме.
7) Рассчитать и построить импульсные характеристики ЦФ. Сравнить
построенные характеристики с аналогичными характеристиками АФпрототипа. Построить логарифмические характеристики ЦФ в функции
абсолютной псевдочастоты.
4
2. Проектирование КИХ – фильтра
1) В соответствии с параметрами ЦФ (таблица 1) провести синтез КИХ –
фильтра с линейной ФЧХ с использованием окна Хэмминга. Синтез ЦФ
провести без использования специализированных функций MatLab синтеза
ЦФ и с их использованием. Этапы синтеза ЦФ отразить в работе.
2) Записать передаточную функцию ЦФ. Построить диаграмму нулей.
3) С помощью системы MatLab построить импульсную и переходную
характеристики фильтра и сигнал на выходе для произвольного входа
(входной сигнал необходимо задавать так, чтобы показать правильность
проведенного синтеза ЦФ).
4) Рассчитать и построить частотные характеристики ЦФ.
5) Составить структурную схему фильтра в последовательной или
параллельной форме.
3. Проектирование и анализ фильтров в графической среде FDATool
Используя графическую среду проектирования FDATool, рассчитать и
получить характеристики проектируемых дискретных фильтров. Сравнить
полученные результаты.
5
Проектирование КИХ – фильтра
Поскольку отсчеты импульсной характеристики КИХ – фильтра
одновременно являются и коэффициентами его передаточной функции,
задача синтеза сводится к получению импульсной характеристики фильтра.
Требования к фильтра приведены в таблице 1.
Вычислим импульсную характеристику hи(n) идеально фильтра.
Обозначим частоту разрыва идеальной АЧХ через 0 . Тогда для 0
получим неравенство  p  0  s , при ширине переходной полосы
  s   p . Удобно определить 0 такой, чтобы  p и s располагались
симметрично относительно 0 . Иначе говоря, будем полагать, что
  s
.
0  p
2
В случае ПФ, расчет должен выполняться с учетом двух переходных
полос:
  minp1  s1 , s1  p1 .
s1  2    Fs1  2    3000  6000 ;
p1  2    Fp1  2    5000  10000 ;
s 2  2    Fs 2  2   11000  22000 .
p2  2    Fp2  2    9000  18000 ;
  min p1  s1  , s 2  p 2   min 10000  6000  ,  22000  18000   4000 .
Выражение для импульсной характеристики избирательного фильтра
при усечении до N членов (причем N – нечетное):
 T sin02nT d 01T d sin01nTd
 
h(n)  02 d

;
01  p1 s1  8000 ;

02nTd

01nTd
2
02Td 01T d
  s 2
h(0) 

;
02  p 2
 20000 .


2
20000 n
8000 n
sin
20000
50000  8000
50000  0, 4 sin 0, 4 n  0,16 sin 0,16 n 
h(n) 
50000 20000 n
50000 8000 n
0, 4 n
0,16 n
50000
50000
sin 0, 4 n  sin 0,16 n

;
n
sin
h(0)  0,4  0,16  0,24.
Выберем порядок фильтра R=N-1.
Ширина переходной полосы, вследствие свертки в частотной области
прямоугольной характеристики и окна равно ширине главного лепестка окна,
величина которого может быть определена для некоторого типов окон. Это
означает, что, во-первых, чем более узкую переходную полосу требуется
получить, тем больше должна быть длина окна N, а потому и порядок
фильтра; во-вторых, чем большее подавление требуется в полосе
6
задерживания, тем более гладкое окно необходимо использовать; последнее
приводит к увеличению ширины переходной полосы и, как следствие, к
увеличению порядка фильтра с целью выполнения заданных требований. В
связи с этим при проектировании ЦФ следует придерживаться
рекомендаций: во-первых, подобрать окно с подходящим уровнем
пульсаций; во-вторых, выбрать число отсчетов N, обеспечивающее
требуемую ширину переходной полосы.
Если окно выбрано, то N для окон Хэмминга оценивается как
N  k T  ,
д 
где k=4;  - ширина переходной полосы.
4
N
1
4000
50000
 50.
Так как, N должно быть нечетным, примем N=49. Тогда порядок
фильтра R=N-1=48.
Порядок фильтра слишком высок, изменим исходные параметры.
Примем:
Таблица 2
АФ-прототип
ЦФ
Мето
Частота
д
Тип
дискретизаци
Rp ,
Rs ,
Fp , кГц
Fs , кГц
синте
фильтр АЧХ
и f д , кГц
дБ
дБ
за
а
ПФ
Б
1,5
40
ИП
50
7,0…8,0 2,0…13,0
p1  2    Fp1  2    7000  14000 ;
p2  2    Fp2  2    8000  16000 ;
s1  2    Fs1  2    2000  4000 ;
s 2  2    Fs 2  2   13000  26000 .
  min p1  s1  , s 2  p 2   min 14000  4000  ,  26000  16000   10000 .
N
4
1
10000
50000
 20.
Тогда, порядок фильтра R=N-1=19.
Вычислим N отсчетов функции окна Хэмминга:
w n  0,54  0,46cos
7
2 n
,
N 1
Отсчеты весовой функции окна Хэмминга
Рассчитаем импульсную характеристику реального фильтра по формуле
h(n)  hи (n)w(n).
8
Импульсная характеристика идеального фильтра
Проверим выполнение заданных требований
Рассчитаем АЧХ и ЛАЧХ
9
Из ЛАЧХ видно, что полосы пропускания и задержки не соответствуют
заданным требованиям.
Увеличение порядка фильтра и изменение ширины переходной полосы
не дает положительных результатов. В соответствии с указаниями к курсовой
работе выберем другое окно и повторим процедуру.
Для обеспечения заданного минимального затухания в полосе
задерживания и малой пульсации в полосе пропускания необходимо
подобрать окно с подходящим уровнем пульсаций и выбрать длину N окна,
при которой обеспечивается требуемая переходная полоса пропускания.
Однако это приводит к излишнему росту порядка фильтра и трудностями с
его реализацией. Эта проблема решается с помощью окна Кайзера.
Окно Кайзера определяется формулой:
I ( )
w(n)  0
, при 0  n  N 1 ,
I 0 ( )
где  - параметр окна, определяющий величину пульсации
2n  N  1 

 N 1 
  1  
2
I 0 ( x) - функция Бесселя первого рода нулевого порядка
Определим   min(1 , 2 ) , где
10
0,05 R

10 p 1 100,051,5
1
1  0,05Rp  0,051,5
 0,085 ,

1
10
1 10
2  100,05R  100,0540  0,01.
S
Тогда,   0,01
Для найденного  , определим Rs
Rs  20lg( )  20lg(0,01)  40
Определим параметр  на основании эмпирических выражений
0,

  0,5842( Rs  21)0,4  0,07886( Rs  21),
0,1102( R  8,7),
s

и параметр D
0,922,
D
(Rs  7,95) / 14,36,
  3,395 , D = 2,232 .
Rs  21;
21  Rs  50;
Rs  50.
Rs  21;
Rs  21.
Выберем наименьшее нечетное значение числа
удовлетворяющее неравенству:
D
2    50000  2,232
N  д 1 
 1  23 .

10000  
Отсчеты весовой функции окна Кайзера
11
отсчетов
N,
Импульсная характеристика идеального фильтра
Импульсная характеристика реального фильтра
12
Для проверки заданных требований построим АЧХ и ЛАЧХ.
Для проверки был произведен расчет фильтра MatLab-функцией fir1, как
видно из АЧХ результаты расчетов одинаковы.
b = fir1(N,wnp,'',wk,'noscale');
N – порядок рассчитываемого фильтра
wnp – двухэлементный вектор частот среза
wk – отсчеты весовой функции
13
Увеличенное ЛАЧХ в области полосы пропускания
Из ЛАЧХ видно что полосы пропускания и задержки удовлетворяют
заданным требованиям.
Запишем передаточную функцию ЦФ по импульсной характеристике и
построим диаграмму нулей.
H ( z)  -0,00119 + 0,00897z1 + 0,00807z2 + 0,00031z3 + 0,00955z4 + 0,03795z5 + 0,03277z6 - 0,05093z7 
- 0,15221z8 - 0,13614z9 + 0,03472z10 + 0,212z11 + 0,212z12 + 0,03472z13 - 0,13614z14 - 0,15221z15 - 0,0509z16 + 0,03277z17 + 0,03795z18 + 0,00955z19 + 0,0003z20 + 0,00807z21 + 0,00897z22 + 0,00119z23
14
Диаграмма нулей
Импульсная и переходная характеристики фильтра
impz(hk,1); - импульсная характеристика
p = filter(hk,1,[0 ones(1, length(hk)-1)]); - переходная характеристика
hk – импульсная характеристика реального фильтра
15
Пропустим через фильтра тестовые сигналы частотами лежащими в
полосе задерживания, в полосе пропускания и в переходной полосе.
Подав на вход фильтра с тестовый сигнал с единичной амплитудой и с
частотой 1500 Гц, видим, что входной сигнал подавлен.
Подав на вход сигнал с частотой 7500 Гц, видим, что входной сигнал не подавлен, т.
к. находится в полосе пропускания, но смещен по фазе на 180 градусов.
16
Подав на вход фильтра сигнал с частотой 14000 Гц, видим, что входной сигнал
подавлен.
Из графиков видно, что в полосе задерживания сигнал практически
полностью подавляется.
Сигнал, с частотой, лежащей в переходной полосе имеет несколько
меньшую амплитуду.
17
Частотные характеристики фильтра
[H, W] = freqz(hk,[1]);
plot(W*Fd/(2*pi),20*log10(abs(H)));
faza = angle(H); faza = unwrap(faza);
plot(W*Fd/(2*pi),faza);
H – вектор комплексной АЧХ;
W – вектор частот;
Функция P = angle(Z) возвращает массив значений аргументов для элементов Z;
Функция P = unwrap(P) корректирует фазовые углы элементов, чтобы убрать
разрывы функции.
Составим структурную схему фильтра в параллельной форме.
Y ( z)
H ( z) 
 -0,00119 + 0,00897z1 + 0,00807z2 + 0,00031z3 + 0,00955z4 + 0,03795z5 + 0,03277z6 - 0,05093z7 
X ( z)
- 0,15221z8 - 0,13614z9 + 0,03472z10 + 0,212z11 + 0,212z12 + 0,03472z13 - 0,13614z14 - 0,15221z15 - 0,0509z16 + 0,03277z17 + 0,03795z18 + 0,00955z19 + 0,0003z20 + 0,00807z21 + 0,00897z22 + 0,00119z23
Y ( z)  X ( z)  H ( z)
X ( z)  xk ,
Y ( z)  yk ,
z n  xk  xk n .
yk  -0,00119xk + 0,00897xk 1 + 0,00807xk 2 + 0,00031xk 3 + 0,00955xk 4 + 0,03795xk 5 +
+ 0,03277xk 6 - 0,05093xk 7 - 0,15221xk 8 - 0,13614xk 9 + 0,03472xk 10 + 0,212xk 11 +
+ 0,212xk 12 + 0,03472xk 13 - 0,13614xk 14 - 0,15221xk 15 - 0,0509xk 16 + 0,03277xk 17 +
0,03795xk 18 + 0,00955xk 19 + 0,0003xk 20 + 0,00807xk 21 + 0,00897xk 22 + 0,00119xk 23
Коэффициенты перед xk , xk 1 ,...,xk 19 обозначим через b0,b1,…,b23 соответственно
18
Тогда структурная схема фильтра в параллельной форме будет иметь вид
изображенный на рисунке.
b0
xk

b1
z 1
z
1
b2
…
b23
z 1
19
yk
Проектирование БИХ – фильтра
При синтезе частотно-избирательных рекурсивных фильтров (ФНЧ,
ФВЧ, ПФ, РФ) удобнее всего воспользоваться хорошо развитым аппаратом
расчета АФ и методами отображения p-плоскости в z-область, т.е. методами
преобразования АФ в ЦФ.
Такой синтез включает в себя:
– выбор метода отображения p-плоскости в z-область;
– расчет АФ по требованиям, заданным к ЦФ;
– применение к АФ выбранного метода отображения p-плоскости в zобласть.
Рассчитываемый по требованиям, заданным к ЦФ, АФ называется
фильтром-прототипом, или просто прототипом.
Основными ограничениями для синтеза ЦФ по прототипу являются:
– сохранение существенных АЧХ прототипа в АЧХ ЦФ, что означает
необходимость отображения мнимой оси i p-плоскости на единичную
окружность z-области;
– устойчивый прототип должен быть преобразован в устойчивый ЦФ,
что означает необходимость отображения полюсов устойчивого прототипа из
левой p-полуплоскости внутрь единичного круга z-области.
БИХ - фильтр рассчитывается на основе аналоговых ФНЧ - прототипов
по следующей методике:
1. На основании требований к ЦФ формулируются требования к
соответствующему аналоговому фильтру-прототипу.
2. Для аналогового фильтра-прототипа определяется нормированный
аналоговый ФНЧ-прототип с частотой среза (границей полосы пропускания)
0  1рад/ с и рассчитывается его нормированные полюсы и нули, исходя из
квадрата его АЧХ.
3. С помощью формул преобразования частот производится
преобразование нормированного ФНЧ-прототипа в аналоговый прототип,
соответствующий исходному ЦФ.
4. Производится пересчет нулей и полюсов из аналоговой области в
цифровую.
Идеальная (прямоугольная) форма АЧХ не может быть физически
реализована, поэтому в теории аналоговых фильтров разработан ряд методов
аппроксимации прямоугольных АЧХ. Аппроксимацией фильтра будем
называть передаточную функцию, у которой АЧХ A() приближается к
одной из идеальных характеристик. Такая передаточная функция должна
удовлетворять следующим условиям:
1) она должна быть рациональной функцией от p с вещественными
коэффициентами;
2) ее полюсы (нули знаменателя) должны лежать в левой полуплоскости
p-плоскости;
20
3) степень полинома числителя должна быть равной или меньше степени
полинома знаменателя.
Согласно методике расчета БИХ-фильтра, требования, предъявляемые у
ЦФ, перекладываются на соответствующий аналоговый ФНЧ-прототип.
Рассчитав ФНЧ, можно несложными преобразованиями изменить его частоту
среза, превратить его в ФВЧ, полосовой, либо режекторный фильтр с
заданными параметрами. Вследствие чего расчет аналогового фильтра
начинается с расчета так называемого фильтра-прототипа, представляющего
собой ФНЧ с частотой среза, равной 1 рад/с.
Рассчитаем АФ-прототип с передаточной функцией G( p)
требованиям, заданным к ЦФ, используя функции buttord и butter.
по
[N W0] = buttord(Wp,Ws,Rp,Rs,'s');
Wp,Ws – двухэлементные вектора границы полос пропускания и задерживания;
Rp - допустимый уровень пульсаций АЧХ в полосе пропускания;
Rs - минимально необходимое затухание в полосе задерживания;
's' – признак аналогового расчета;
N – минимально необходимый порядок фильтра;
W0 – частота среза фильтра.
[b, a] = butter(N,W0,'bandpass','s');
'bandpass' – тип фильтра – полосовой фильтр;
b, a – коэффициенты передаточной функции аналогового фильтра.
b = [0 0 0 1,6292e  12 0 0 0],
a  [1 2,35e+4 6,91e+9 1,05e+14 1,52e+19 1,15e+23 1,08e+28]
Передаточная функция АФ-прототипа
G( p) 
1,6292e  12
p + 2,35e+4p + 6,91e+9p + 1,05e+14p3 + 1,52e+19p2 + 1,15e+23p + 1,08e+28
6
5
4
.
21
Частотные характеристики АФ-прототипа
22
Нули и полюсы АФ-прототипа
23
Представим передаточную функцию G( p) в виде суммы простейших
дробей G j ( p) 
Aj
, используя функцию residue.
p  pj
[r, p, k] = residue(b,a);
b,a – коэффициенты передаточной функции;
r вектор-столбец вычетов;
p вектор-столбец полюсов;
k - вектор-строка целой части дробно-рациональной функции.
3155,8  2063,6 j
3155,8  2063,6 j
5883, 4  742 j


p  (3259  52299 j) p  (3259  52299 j) p  (5883  46650 j)
5883, 4  742 j
2727,5  1333, 2 j
2727,5  1333, 2 j



p  (5883  46650 j) p  (2624  42109 j) p  (2624  42109 j)
Определим ИХ g j (t ) каждого звена, соответствующего G j ( p) :
G( p) 
g j (t )  Aj e j .
pt
g1 (t)  (3155  2063,6 j)e(325992299 j )t ; g2 (t)  (3155,8  2063,6 j)e(325952299 j )t ;
g3 (t)  (5883,4  742 j)e(588346650 j )t ; g4 (t)  (5883,4  742 j)e(588346650 j )t ;
g5 (t)  (2727,5 1333,2 j)e(262442109 j )t ; g6 (t)  (2727,5 1333,2 j)e(262442109 j )t .
Проведем их дискретизацию, заменой t  nTd , получив при этом
импульсные характеристики hi (n) звеньев ЦФ.
g1 (nTd )  (3155  2063,6 j)e(325992299 j )nTd ;
g2 (nTd )  (3155,8  2063,6 j)e(325952299 j )nTd ;
g3 (nTd )  (5883,4  742 j)e(588346650 j )nTd ;
g4 (nTd )  (5883,4  742 j)e(588346650 j )nTd ;
g5 (nTd )  (2727,5  1333,2 j)e(262442109 j )nTd ;
g6 (nTd )  (2727,5 1333,2 j)e(262442109 j )nTd .
Передаточную функцию ЦФ определим как сумму отдельно взятых
N
передаточных функций звеньев, т.е. H ( z)   H j ( z) , где H j ( z) 
j 1
Aj
1 e
p jTd
z 1
.
3155  2063,6 j
3155,8  2063,6 j
2727,5  1333,2 j

 ... 

( 325992299 j )Td 1
( 325952299 j )Td 1
1 e
z
1 e
z
1  e(262442109 j )Td z 1
(3155  2063,6 j) z (3155,8  2063,6 j) z
(2727,5  1333,2 j) z


 ... 

( 325992299 j )Td
( 325952299 j )Td
z e
z e
z  e(262442109 j )Td
b0  b1 z 1  b1 z 2  b1 z 3  b1 z 4  b1 z 5  b1 z 6

a0  a1 z 1  a1 z 2  a1 z 3  a1 z 4  a1 z 5  a1 z 6
H ( z) 
Для вычисления коэффициентов полиномов воспользуемся функцией
p = ploy(A), где p – коэффициенты характеристического полином, A – его корни.
az = poly(exp(p*1/Fd));
az – коэффициенты знаменателя передаточной функции;
24
p – вектор-столбец полюсов;
1/Fd – период дискретизации.
Коэффициенты числителя в MatLab вычислим следующим образом:
for m = 1:length(p)
if m == 1
p1 = [p(2:length(p))];
elseif m == length(p)
p1 = [p(1:length(p)-1)];
else
p1 = [p(1:m-1);p(m+1:length(p))];
end;
A(m,:) = poly(exp(p1*1/Fd));
R(m,:) = r(m)*A(m,:)/Fd;
end;
bz1 = 0;
for m = 1:length(p)
bz1 = bz1 + R(m,:);
end;
bz = real([bz1,0]);
Результат этих действий (коэффициенты числителя и знаменателя
передаточной функции ЦФ):
b = [-6,93e-17 0,0044 -0,0093 0,0018 0,0070 -0,0038 0]
a = [1 -3,261 6,086 -6,833 5,202 -2,381 0,624]
Проверка расчета функцией MatLab impinvar дает примерно равный
результат:
[bzp,azp] = impinvar(b,a,Fs);
b,a – коэффиценты числителя и знаменателя функции передачи аналогового
прототипа;
Fs – частота дискретизации;
bzp,azp – коэффициенты числителя и знаменателя функции передачи ЦФ.
b = [-6,36e-17 0,0043
0,0093 0,0017 0,0069 -0,0037 0]
a = [1 -3,2611 6,0865 -6,8335 5,2023 -2,3813 0,6246]
H ( z) 
6,36 1017  0,0043z1  0,0093z 2  0,0017 z 3  0,0069z 4  0,0037 z 5
1  3,2611z 1  6,0865z 2  6,8335z 3  5,2023z 4  2,3813z 5  0,6246z 6
25
Построим диаграмму полюсов и нулей
Увеличенная диаграмма нулей в области большего числа нулей
Все полюсы фильтра лежат внутри окружности единичного радиуса zплоскости, следовательно, фильтр устойчив.
Составим разностное уравнение ЦФ и решим его методом итерационной
процедуры для управляющего сигнала типа единичный скачок.
26
y(k )  3,2611y(k  1)  6,0865 y(k  2)  6,8335 y(k  3)  5,2023 y(k  4) 
2,3813 y(k  5)  0,6246 y(k  6)  6,36*1017 x(k )  0,0043x(k 1) 
0,0093x(k  2)  0,0017 x(k  3)  2,3813x(k  4)  0,6246 x(k  5)
Решение методом итерационной процедуры в MatLab:
for n = 1:length(h1)
y(n) = 0;
for m=1:n
y(n)=y(n)+x(m)*h1(n-m+1);
end
end;
где
h1 – импульсная характеристика ЦФ;
x – вектор управляющего сигнала типа единичный скачок.
Решение разностного уравнения для управляющего сигнала типа
единичный скачок
27
Импульсная и переходная характеристики фильтра
Построим сигналы на выходе фильтра, подавая на вход сигналы с
частотой лежащей в полосе задержки и в полосе пропускания.
Подав на вход фильтра с тестовый сигнал с единичной амплитудой и с частотой 1500
Гц, видим, что входной сигнал подавлен.
28
Подав на вход сигнал с частотой 7500 Гц, видим, что входной сигнал не подавлен, т.
к. находится в полосе пропускания.
Подав на вход фильтра сигнал с частотой 14000 Гц, видим, что входной сигнал
подавлен.
29
Из графиков видно, что в полосе задерживания сигнал практически
полностью подавляется.
Составим каноническую структурную схемы фильтра. Для этого
запишем уравнения:
17
1
2
3
4
5
Y ( z)  6,36 10  0,0043z  0,0093z  0,0017 z  0,0069z  0,0037 z U ( z)
H ( z) 

;
X ( z) 1  3,2611z 1  6,0865z 2  6,8335z 3  5,2023z 4  2,3813z 5  0,6246z 6 U ( z)
Y ( z)   6,36 1017  0,0043z 1  0,0093z 2  0,0017 z 3  0,0069z 4  0,0037 z 5 U ( z);
U ( z)  X ( z)   3,2611z 1  6,0865z 2  6,8335z 3  5,2023z 4  2,3813z 5  0,6246z 6 U ( z);
yk  6,36 1017 uk  0,0043uk 1  0,0093uk 2  0,0017uk 3  0,0069uk 4  0,0037uk 5 ;
xk  uk  3,2611uk 1  6,0865uk 2  6,8335uk 3  5,2023uk 4  2,3813uk 5  0,6246uk 6 .
xk
uk
b0
a1
z 1
b1
a2
z 1
b2
a3
z 1
b3
a4
z 1
b4
a5
z 1
b5

a6
z 1
30

yk
где b0= 6,36 1017 ; b1= 0,0043 … b5= 0,0037 ;
a0=1 ; a1= 3,2611 … a6= 0,6246 .
Составим параллельную структурную схему фильтра. Для этого
передаточную функцию представим в виде суммы простых дробей с
помощью matlab-функции residue.
[r, p, k] = residue(bz,az);
bz,az – коэффициенты передаточной функции;
r вектор-столбец вычетов;
p вектор-столбец полюсов;
k - вектор-строка целой части дробно-рациональной функции.
H ( z) 
r1
r2
rn


...

 k0  k1z 1  ...  km z (mn)
1
1
1
1  p1z
1  p2 z
1  pn z
-0,0533 
-0,0533 


0,00383
r
;
0,00383


0,05168


0,05168
0,63179 
0,63179 


0,46943 
p
 ;
0,46943


0,52933 


0,52933 
31
k  -6,36e-17 .
xm
r1

p1
r2

p2
r6
z 1
z 1

z 1
p6
k
32

ym
Проектирование и анализ фильтров в графической среде FDATool
Рабочее окно в FDATool (БИХ-фильтр)
АЧХ и ФЧХ БИХ-фильтра
33
Импульсная характеристика
34
Переходная характеристика
Все характеристики, построенные в FDATool для БИХ-фильтра совпадают с
рассчитанными ранее.
35
Рабочее окно в FDATool (КИХ-фильтр)
АЧХ и ФЧХ КИХ-фильтра
36
Импульсная характеристика
Переходная характеристика
37
Диаграмма нулей
Все характеристики, построенные в FDATool для КИХ-фильтра совпадают с
рассчитанными ранее.
38
Листинг MatLab-программы
KIH.m
%Проектирование КИХ фильтра (ПФ)
clc;
clear all;
close all;
global Fn;
% Исходные данные
% Синтез с использованием окна Хэмминга
% Допустимый уровень пульсаций АЧХ в полосе пропускания
Rp = 1.5; %[ДБ]
% Минимально необходимое затухание АЧХ в полосе задерживания
Rs = 40;
%[ДБ]
% Границы полосы пропускания
Fp1 = 7000;
Fp2 = 8000; % [Гц]
Wp1 = 2*pi*Fp1;
% [рад/c]
Wp2 = 2*pi*Fp2;
% [рад/c]
% Границы полосы задерживания
Fs1 = 2000; % [Гц]
Fs2 = 13000; % [Гц]
Ws1 = 2*pi*Fs1;
% [рад/c]
Ws2 = 2*pi*Fs2;
% [рад/c]
Fd = 50000; % Частота дискретизации % [Гц]
Td = 1/Fd;
Wd = 2*pi*Fd;
Fn = Fd/2; % Частота Найквиста
%Частоты среза [рад/c]
W01 = (Wp1 + Ws1)/2;
W02 = (Wp2 + Ws2)/2;
%Частоты для тестирования фильтра
f1 = 1500;
f2 = 7500;
f3 = 14000;
%--- Окно Хэмминга
dW = min((Wp1-Ws1),(Ws2-Wp2));
N = round(4*pi*Fd/dW);
if mod(N,2)==0
N = N+1;
end;
n = 0:N-1;
R = N-1;
%-- выражение для импульсной характериски идеального фильтра
i = 0;
hn(N)=(W02-W01)*Td/pi;
39
for k = -(N-1)/2:(N-1)/2;
i = i+1;
if k~=0
hn(i) = (W02*Td/pi)*(sin(W02*k*Td)/(W02*k*Td))(W01*Td/pi)*(sin(W01*k*Td)/(W01*k*Td));
else
hn(i) = (W02*Td/pi)-(W01*Td/pi);
end;
end;
figure; stem(hn); grid;
title('Импульсная характеристика идеального фильтра окна
Хэмминга ');
%-- отсчеты весовой функции w(n)
w = 0.54-0.46*cos(2*pi*n/R);
figure; stem(w,'bo'); grid on; hold;
xlabel('n');
ylabel('w(n)');
title('Отсчеты весовой функции окна Хэмминга');
% Вычисление отсчетов функции окна средствами MatLab
wn = hamming(N);
stem(wn,'b.'); grid on;
xlabel('n');
ylabel('w(n)');
legend('полученные по формуле окна Хэмминга','полученные
функцией hamming(N)');
% Вычисление импульсной характеристики реального фильтра
h = hn.*wn';
figure; stem(h); grid;
title('Импульсная характеристика реального фильтра окна
Хэмминга');
% Частотные характеристики
[H W] = freqz(h,[1]);
figure; plot(W*Fd/(2*pi),abs(H),'b'); grid; hold on;
title('АЧХ (окно Хэмминга)');
xlabel('f, Hz');
ylabel('A');
figure; plot(W*Fd/(2*pi),20*log10(abs(H)), 'b'); grid on; hold
on;
title('ЛАЧХ (окно Хэмминга)');
xlabel('f, Hz');
ylabel('A, dB');
%-- построим границы полос пропускания и задержки
plot([Fs1, Fs1], [0, -Rs], 'r');
plot([Fs2, Fs2], [0, -Rs], 'r');
plot([0, Fs1] ,[-Rs, -Rs], 'r');
plot([Fs2, Fs2+1000] ,[-Rs, -Rs], 'r');
plot([Fp1, Fp1], [0, -Rp], 'r');
plot([Fp2, Fp2], [0, -Rp], 'r');
plot([Fp1, Fp2], [-Rp, -Rp], 'r');
40
%--- Окно Кайзера
%-- определим delta_min
d1 = (10^(0.05*Rp)-1)/(10^(0.05*Rp)+1);
d2 = 10^(-0.05*Rs);
d = min(d1, d2);
Rsn = -20*log10(d);
if Rsn<=21;
alpha = 0;
elseif Rsn<=50;
alpha = 0.5842*(Rsn-21)^0.4 + 0.07886*(Rsn-21);
elseif Rsn>50;
alpha = 0.1102*(Rsn-8.7);
end;
if Rsn<=21;
D = 0.992;
else
D = (Rsn-7.95)/14.36;
end;
Nk = ceil(Wd*D/dW + 1);
if mod(Nk,2)==0
Nk = Nk+1;
end;
Nk = Nk - 1;
%-- импульсная харктеристика идеального фильтра
i = 0;
hi(Nk)=(W02-W01)*Td/pi;
for k = -(Nk-1)/2:(Nk-1)/2;
i = i+1;
if k~=0
hi(i) = (W02*Td/pi)*(sin(W02*k*Td)/(W02*k*Td))(W01*Td/pi)*(sin(W01*k*Td)/(W01*k*Td));
else
hi(i) = (W02*Td/pi)-(W01*Td/pi);
end;
end;
figure; stem(hi); grid;
title('Импульсная характеристика идеального фильтра окна
Кайзера');
%-- отсчеты весовой функции
wk = kaiser(Nk,alpha)';
figure; stem(wk,'b.'); grid;
xlabel('n');
ylabel('w(n)');
title('Отсчеты весовой функции окна Кайзера');
%-- ипульсная характеристика реального фильтра
hk = hi.*wk;
figure; stem(hk); grid;
title('Импульсная характеристика реального фильтра окна
Кайзера');
%Проверка выполнения заданных требований
%АЧХ, ЛАЧХ
[Hk Wk] = freqz(hk,[1]);
41
figure; plot(Wk*Fd/(2*pi),abs(Hk),'b'); grid on; hold;
title('АЧХ (окно Кайзера)');
xlabel('f, Hz');
ylabel('A');
%Синтез с использованием специализированных функций MatLab
b = fir1(N-1, [W01/(2*pi*Fd/2) W02/(2*pi*Fd/2)], hamming(N));
[Hfir Wfir] = freqz(b,[1]);
plot(Wfir*Fd/(2*pi),abs(Hfir),'m:');
legend('АЧХ фильтра без спец. функций','АЧХ фильтра Matlab
функцией');
%
figure; plot(Wk*Fd/(2*pi),20*log10(abs(Hk)), 'b'); grid on; hold
on;
title('ЛАЧХ (окно Кайзера)');
xlabel('f, Hz');
ylabel('A, dB');
%построим границы полос пропускания и задержки
plot([Fs1, Fs1], [0, -Rsn], 'r');
plot([Fs2, Fs2], [0, -Rsn], 'r');
plot([0, Fs1] ,[-Rsn, -Rsn], 'r');
plot([Fs2, Fs2+1000] ,[-Rsn, -Rsn], 'r');
plot([Fp1, Fp1], [0, -Rp], 'r');
plot([Fp2, Fp2], [0, -Rp], 'r');
plot([Fp1, Fp2], [-Rp, -Rp], 'r');
%построим диаграмму нулей
figure; zplane(hk,[1]); grid;
title ('Диаграмма полюсов и нулей');
%импульсная и переходная характеристики фильтра
figure; subplot(2,1,1); impz(hk,1); grid;
title('Импульсная характеристика');
p = filter(hk,1,[0 ones(1, length(hk)-1)]);
subplot(2,1,2); stem((0:length(p)-1), p); grid;
title('Переходная характеристика');
% Прохождение тестового сигнала через фильтр
[y1, x1, t1] = testfilter (h, 1, Fd, f1, 200);
[y2, x2, t2] = testfilter (h, 1, Fd, f2, 100);
[y3, x3, t3] = testfilter (h, 1, Fd, f3, 50);
% Прохождение тестовых сигналов через фильтр
figure;
plotxy (t1*1e3, y1, 'r', 't, мc', 'y', 'Тестирование фильтра для
1500 Гц', 'plot'); hold on;
plot (t1*1e3, x1, 'k--');
legend ('Выходной сигнал для 1500 Гц','Исходный сигнал'), hold
off;
figure;
plotxy (t2*1e3, y2, 'r', 't, мc', 'y', 'Тестирование фильтра для
7500 Гц', 'plot'); hold on;
plot (t2*1e3, x2, 'k--');
42
legend ('Выходной сигнал для 7500 Гц','Исходный сигнал'), hold
off;
figure;
plotxy (t3*1e3, y3, 'r', 't, мc', 'y', 'Тестирование фильтра для
14000 Гц', 'plot'); hold on;
plot (t3*1e3, x3, 'k--');
legend ('Выходной сигнал для 14000 Гц','Исходный сигнал'), hold
off;
%АЧХ, ФЧХ
[H, W] = freqz(hk,[1]);
figure; subplot(2,1,1); plot(W*Fd/(2*pi),20*log10(abs(H)));
grid;
xlabel('f, Hz'); ylabel('Magnituda'); title('AChH');
faza = angle(H); faza = unwrap(faza);
subplot(2,1,2); plot(W*Fd/(2*pi),faza); grid;
xlabel('f, Hz'); ylabel('faza, grad'); title('FChH');
% Последовательная структурная схема фильтра
[sos, G] = tf2sos (h, 1) %преобразование коэффициентов полиномов
функции передачи в параметры набора секций второго порядка
BIH.m
clc; clear all; close all;
Fs = 50000;
Rp = 1.5; Rs = 40;
Fp1 = 7000; Fp2 = 8000;
Fs1 = 2000; Fs2 = 13000;
Fd = 50000; Td = 1/Fd;
Wd = 2*pi*Fd;
Wp1 = 2*pi*Fp1; Wp2 = 2*pi*Fp2;
Ws1 = 2*pi*Fs1; Ws2 = 2*pi*Fs2;
dW = min((Wp1-Ws1),(Ws2-Wp2));
%Частоты для тестирования фильтра
f1 = 1500;
f2 = 7500;
f3 = 14000;
%--- АФ-прототип
Wp = [Wp1, Wp2];
Ws = [Ws1, Ws2];
[N W0] = buttord(Wp,Ws,Rp,Rs,'s');
[b, a] = butter(N,W0,'bandpass','s');
w = 0:10:Wd/2;
h = freqs(b,a,w);
figure;
subplot(2,1,1); plot(w/(2*pi),abs(h)); grid on; hold;
plot([0 W0(1)/(2*pi) W0(1)/(2*pi) W0(2)/(2*pi) W0(2)/(2*pi)
Wd/(4*pi)], [0 0 1 1 0 0],'r--');
legend('АЧХ АФ-прототипа','АЧХ идеального фильтра');
title('АЧХ');
xlabel('f,Hz'); ylabel('Magnituda');
subplot(2,1,2); plot(w/(2*pi),unwrap(angle(h))*180/pi); grid;
43
title('АЧХ'); xlabel('f, Hz'); ylabel('фаза, grad');
figure; plot(w/(2*pi),20*log10(abs(h))); grid on; hold on;
title('ЛАЧХ'); xlabel('f, Hz'); ylabel('Magnituda, Db');
plot([Fs1, Fs1], [0, -Rs], 'r');
plot([Fs2, Fs2], [0, -Rs], 'r');
plot([0, Fs1] ,[-Rs, -Rs], 'r');
plot([Fs2, Fs2+1000] ,[-Rs, -Rs], 'r');
plot([Fp1, Fp1], [0, -Rp], 'r');
plot([Fp2, Fp2], [0, -Rp], 'r');
plot([Fp1, Fp2], [-Rp, -Rp], 'r');
figure; [hz, hp, hl] = zplane(b,a); grid;
%-- импульсная характеристика АФ-прототипа
haf = impz(b,a,150);
figure; stem(haf); grid; title('Импульсная характеристика АФпрототипа');
%-- представим передаточную функцию G(p) в виде простейших
дробей
[r, p, k] = residue(b,a);
%-- вычислим импульсные характеристики АФ и ЦФ
%-- для каждого звена, задем дискретизируем
Tm = 0.002;
taf = 0:1/Fd*1e-2:Tm;
ncf = 0:1/Fd:Tm;
for j = 1:length(r)
g(j,:) = r(j)*exp(p(j)*taf);
hn(j,:) = r(j)*exp(p(j)*ncf);
end;
%-- передаточная функция ЦФ
az = poly(exp(p*1/Fd));
for m = 1:length(p)
if m == 1
p1 = [p(2:length(p))];
elseif m == length(p)
p1 = [p(1:length(p)-1)];
else
p1 = [p(1:m-1);p(m+1:length(p))];
end;
A(m,:) = poly(exp(p1*1/Fd));
R(m,:) = r(m)*A(m,:)/Fd;
end;
bz1 = 0;
for m = 1:length(p)
bz1 = bz1 + R(m,:);
end;
bz = real([bz1,0]);
%-- проверка средствами Matlab
[bzp,azp] = impinvar(b,a,Fs);
%-- построим диаграмму полюсов и нулей
figure; zplane(bzp,azp); grid;
%-- решим разностное уравнение методом итерационной процедуры
%для управляющего сигнала типа единичный скачок
h1 = impz(bzp,azp,150);
x = [ones(1,1),zeros(1,length(h1)-1)];
44
for n = 1:length(h1)
y(n) = 0;
for m=1:n
y(n)=y(n)+x(m)*h1(n-m+1);
end
end;
figure; plot(0:length(x)-1, y(1:length(x))); grid;
xlabel('n'); ylabel('amplituda');
%-- построим импульсную и переходную характеристики
figure; per = filter(bzp, azp, ones(1, length(h1)));
subplot(2,1,1); stem(h1); grid;
title('Импульсная характеристика');
xlabel('n'); ylabel('amplituda');
subplot(2,1,2); stem(per); grid;
title('Переходная характеристика');
xlabel('n'); ylabel('amplituda');
%-- тестовые сигналы
% Прохождение тестового сигнала через фильтр
[y1, x1, t1] = testfilter (bzp, azp, Fd, f1, 200);
[y2, x2, t2] = testfilter (bzp, azp, Fd, f2, 100);
[y3, x3, t3] = testfilter (bzp, azp, Fd, f3, 50);
% Прохождение тестовых сигналов через фильтр
figure;
plotxy (t1*1e3, y1, 'r', 't, мc', 'y', 'Тестирование фильтра для
1500 Гц', 'plot'); hold on;
plot (t1*1e3, x1, 'k--');
legend ('Выходной сигнал для 1500 Гц','Исходный сигнал'), hold
off;
figure;
plotxy (t2*1e3, y2, 'r', 't, мc', 'y', 'Тестирование фильтра для
7500 Гц', 'plot'); hold on;
plot (t2*1e3, x2, 'k--');
legend ('Выходной сигнал для 7500 Гц','Исходный сигнал'), hold
off;
figure;
plotxy (t3*1e3, y3, 'r', 't, мc', 'y', 'Тестирование фильтра для
14000 Гц', 'plot'); hold on;
plot (t3*1e3, x3, 'k--');
legend ('Выходной сигнал для 14000 Гц','Исходный сигнал'), hold
off;
%АЧХ, ФЧХ
[hcf wcf] = freqz(bzp,azp);
figure; plot(wcf*Fd/(2*pi),20*log10(abs(hcf))); grid; hold on;
xlabel('f, Hz'); ylabel('Amplituda(dB)'); title('ЛАЧХ');
plot(w/(2*pi),20*log10(abs(h)),'m:');
legend('ЦФ','АФ-прототип');
plot([Fs1, Fs1], [0, -Rs], 'r');
plot([Fs2, Fs2], [0, -Rs], 'r');
plot([0, Fs1] ,[-Rs, -Rs], 'r');
plot([Fs2, Fs2+1000] ,[-Rs, -Rs], 'r');
plot([Fp1, Fp1], [0, -Rp], 'r');
plot([Fp2, Fp2], [0, -Rp], 'r');
plot([Fp1, Fp2], [-Rp, -Rp], 'r');
45
figure; plot(wcf*Fd/(2*pi),unwrap(angle(hcf))*180/pi); grid;
hold on;
xlabel('f, Hz'); ylabel('фаза, grad'); title('ФЧХ');
plot(w/(2*pi),unwrap(angle(h))*180/pi,'m');
legend('ЦФ','АФ-прототип');
[r,p,k] = residue(bzp,azp);
plotxy.m
% Функция строит графики с подписями
function plotxy (x, y, color, strX, strY, strTitle, str);
if
(strcmp (str, 'stem')), stem (x, y, color)
elseif (strcmp (str, 'plot')), plot (x, y, color)
else
error ('Ошибка в параметре построения графика!');
end;
grid on;
xlabel (strX);
ylabel (strY);
title (strTitle);
testfilter.m
% Прохождение тестового сигнала через фильтр
function [yt, ya, ta] = testfilter (b, a, fd, f, N);
n = 0 : N - 1;
ta = 0 : 1e-6 : (N - 1) / fd;
ya = sin (f*2*pi*n/fd);
yt = pulstran (ta*fd, [n' filter(b, a, sin (f*2*pi*n/fd))'] ,
'sinc');
ya = pulstran (ta*fd, [n' ya'], 'sinc');
46
Заключение
В данной курсовой работе был произведен расчет КИХ- и БИХ- фильтров.
Расчет фильтров осуществлялся как стандартными средствами Matlab, так и с
использованием
пакета
Signal
Processing,
обеспечивающий
широкие
возможности по созданию программ обработки сигналов для современных
научных и технических приложений. Результаты расчетов идентичны,
следовательно, синтез фильтров выполнен правильно.
47
Список использованной литературы
1. Малинин Г.В., Лазарева Н.М. Синтез цифровых фильтров: Руководство
к курсовой работе. – Чебоксары: Изд-во Чуваш. ун-т., 2006. 68 с.
2. Сергиенко А.Б. Цифровая обработка сигналов. - СПб.: Питер, 2002. 608 с.
3. Лазарева Н.М., Антонов В.И. Базовые алгоритмы цифровой обработки
сигналов: Метод. Указания к лабораторным работам. – Чебоксары: Издво Чуваш. ун-т., 2003, 60с.
48
Документ
Категория
Радиоэлектроника
Просмотров
15
Размер файла
1 096 Кб
Теги
работа, курсовая
1/--страниц
Пожаловаться на содержимое документа