Курсовая работа по теории сигналов «Синтез цифровых фильтров» Вариант 15

ФГБОУ ВПО «Чувашский государственный университет

имени »

Кафедра промышленной электроники

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

по теории сигналов

«Синтез цифровых фильтров»

Вариант 15

Выполнил:

ст. гр. РТЭ 31-09

Проверил:

Чебоксары 2012

Содержание

1.  Аннотация ……………………………………………………………....…3

2.  Задание на курсовую работу ……………………………………………..4

3.  Проектирование КИХ-фильтра …………………………………………..6

4.  Проектирование БИХ-фильтра ………………………………..………...20

5.  Проектирование и анализ фильтров в графической среде FDATool.....33

6.  Листинг MatLab-программы …………………………………………….39

7.  Заключение …………………………………………………………...…..47

8.  Список использованной литературы ………………………....................48

Аннотация

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

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

Таблица 1

АФ-прототип	ЦФ	Метод синтеза	Частота дискретизации , кГц
Тип фильтра	АЧХ	, дБ	, дБ	, кГц	, кГц
ПФ	Б	1,5	40	5,0…9,0	3,0…11,0	ИП	50

Используемые сокращения:

Б – фильтр Баттерворта

– допустимый уровень пульсаций АЧХ в полосе пропускания

– минимально необходимое затухание в полосе задержки

, – границы полос пропускания и задержки

ИП – инвариантная импульсная характеристика

1. Проектирование БИХ – фильтра

В соответствии с таблицей 1 необходимо:

1) С помощью системы MatLab синтезировать АФ-прототип. Записать передаточную функцию АФ. Рассчитать и построить частотные характеристики АФ-прототипа. Сравнить АЧХ АФ с АЧХ идеального фильтра.

2) Синтезировать рекурсивный ЦФ по аналоговому прототипу. Синтез ЦФ провести без использования специализированных функций MatLab синтеза ЦФ и с их использованием. Этапы синтеза отразить в работе.

3) Записать передаточную функцию ЦФ. Построить диаграмму полюсов и нулей.

4) Составить разностное уравнение ЦФ. Решить разностное уравнение методом итерационной процедуры для управляющего сигнала типа единичный скачок. Построить кривую переходного процесса. С помощью системы MatLab построить импульсную и переходную характеристики фильтра и сигнала на выходе для произвольного входа.

5) Определит устойчивость фильтра.

6) Составить каноническую структурную схему фильтра. Разложением передаточной функции ЦФ на произведение/сумму звеньев второго (первого) порядка реализовать структурную схему фильтра в последовательной или параллельной форме.

7) Рассчитать и построить импульсные характеристики ЦФ. Сравнить построенные характеристики с аналогичными характеристиками АФ-прототипа. Построить логарифмические характеристики ЦФ в функции абсолютной псевдочастоты.

2. Проектирование КИХ – фильтра

1) В соответствии с параметрами ЦФ (таблица 1) провести синтез КИХ – фильтра с линейной ФЧХ с использованием окна Хэмминга. Синтез ЦФ провести без использования специализированных функций MatLab синтеза ЦФ и с их использованием. Этапы синтеза ЦФ отразить в работе.

2) Записать передаточную функцию ЦФ. Построить диаграмму нулей.

3) С помощью системы MatLab построить импульсную и переходную характеристики фильтра и сигнал на выходе для произвольного входа (входной сигнал необходимо задавать так, чтобы показать правильность проведенного синтеза ЦФ).

4) Рассчитать и построить частотные характеристики ЦФ.

5) Составить структурную схему фильтра в последовательной или параллельной форме.

3. Проектирование и анализ фильтров в графической среде FDATool

Используя графическую среду проектирования FDATool, рассчитать и получить характеристики проектируемых дискретных фильтров. Сравнить полученные результаты.

Проектирование КИХ – фильтра

Поскольку отсчеты импульсной характеристики КИХ – фильтра одновременно являются и коэффициентами его передаточной функции, задача синтеза сводится к получению импульсной характеристики фильтра.

Требования к фильтра приведены в таблице 1.

Вычислим импульсную характеристику hи(n) идеально фильтра.

Обозначим частоту разрыва идеальной АЧХ через . Тогда для получим неравенство , при ширине переходной полосы . Удобно определить такой, чтобы и располагались симметрично относительно . Иначе говоря, будем полагать, что

.

В случае ПФ, расчет должен выполняться с учетом двух переходных полос:

.

Выражение для импульсной характеристики избирательного фильтра при усечении до N членов (причем N – нечетное):

Выберем порядок фильтра R=N-1.

Ширина переходной полосы, вследствие свертки в частотной области прямоугольной характеристики и окна равно ширине главного лепестка окна, величина которого может быть определена для некоторого типов окон. Это означает, что, во-первых, чем более узкую переходную полосу требуется получить, тем больше должна быть длина окна N, а потому и порядок фильтра; во-вторых, чем большее подавление требуется в полосе задерживания, тем более гладкое окно необходимо использовать; последнее приводит к увеличению ширины переходной полосы и, как следствие, к увеличению порядка фильтра с целью выполнения заданных требований. В связи с этим при проектировании ЦФ следует придерживаться рекомендаций: во-первых, подобрать окно с подходящим уровнем пульсаций; во-вторых, выбрать число отсчетов N, обеспечивающее требуемую ширину переходной полосы.

Если окно выбрано, то N для окон Хэмминга оценивается как

,

где k=4; - ширина переходной полосы.

Так как, N должно быть нечетным, примем N=49. Тогда порядок фильтра R=N-1=48.

Порядок фильтра слишком высок, изменим исходные параметры. Примем:

Таблица 2

АФ-прототип	ЦФ	Метод синтеза	Частота дискретизации , кГц
Тип фильтра	АЧХ	, дБ	, дБ	, кГц	, кГц
ПФ	Б	1,5	40	7,0…8,0	2,0…13,0	ИП	50

Тогда, порядок фильтра R=N-1=19.

Вычислим N отсчетов функции окна Хэмминга:

Отсчеты весовой функции окна Хэмминга

Рассчитаем импульсную характеристику реального фильтра по формуле

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

Проверим выполнение заданных требований

Рассчитаем АЧХ и ЛАЧХ

Из ЛАЧХ видно, что полосы пропускания и задержки не соответствуют заданным требованиям.

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

Для обеспечения заданного минимального затухания в полосе задерживания и малой пульсации в полосе пропускания необходимо подобрать окно с подходящим уровнем пульсаций и выбрать длину N окна, при которой обеспечивается требуемая переходная полоса пропускания. Однако это приводит к излишнему росту порядка фильтра и трудностями с его реализацией. Эта проблема решается с помощью окна Кайзера.

Окно Кайзера определяется формулой:

при ,

где - параметр окна, определяющий величину пульсации

- функция Бесселя первого рода нулевого порядка

Определим , где

,

.

Тогда,

Для найденного , определим

Определим параметр на основании эмпирических выражений

и параметр D

, D = 2,232 .

Выберем наименьшее нечетное значение числа отсчетов N, удовлетворяющее неравенству:

.

Отсчеты весовой функции окна Кайзера

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

Импульсная характеристика реального фильтра

Для проверки заданных требований построим АЧХ и ЛАЧХ.

Для проверки был произведен расчет фильтра MatLab-функцией fir1, как видно из АЧХ результаты расчетов одинаковы.

b = fir1(N, wnp,'',wk,'noscale');

N – порядок рассчитываемого фильтра

wnp – двухэлементный вектор частот среза

wk – отсчеты весовой функции

Увеличенное ЛАЧХ в области полосы пропускания

Из ЛАЧХ видно что полосы пропускания и задержки удовлетворяют заданным требованиям.

Запишем передаточную функцию ЦФ по импульсной характеристике и построим диаграмму нулей.

Диаграмма нулей

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

impz(hk,1); - импульсная характеристика

p = filter(hk,1,[0 ones(1, length(hk)-1)]); - переходная характеристика

hk – импульсная характеристика реального фильтра

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

Подав на вход фильтра с тестовый сигнал с единичной амплитудой и с частотой 1500 Гц, видим, что входной сигнал подавлен.

Подав на вход сигнал с частотой 7500 Гц, видим, что входной сигнал не подавлен, т. к. находится в полосе пропускания, но смещен по фазе на 180 градусов.

Подав на вход фильтра сигнал с частотой 14000 Гц, видим, что входной сигнал подавлен.

Из графиков видно, что в полосе задерживания сигнал практически полностью подавляется.

Сигнал, с частотой, лежащей в переходной полосе имеет несколько меньшую амплитуду.

Частотные характеристики фильтра

[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) корректирует фазовые углы элементов, чтобы убрать разрывы функции.

Составим структурную схему фильтра в параллельной форме.

, , .

Коэффициенты перед обозначим через b0,b1,…,b23 соответственно

Тогда структурная схема фильтра в параллельной форме будет иметь вид изображенный на рисунке.

Проектирование БИХ – фильтра

При синтезе частотно-избирательных рекурсивных фильтров (ФНЧ, ФВЧ, ПФ, РФ) удобнее всего воспользоваться хорошо развитым аппаратом расчета АФ и методами отображения p-плоскости в z-область, т. е. методами преобразования АФ в ЦФ.

Такой синтез включает в себя:

– выбор метода отображения p-плоскости в z-область;

– расчет АФ по требованиям, заданным к ЦФ;

– применение к АФ выбранного метода отображения p-плоскости в z-область.

Рассчитываемый по требованиям, заданным к ЦФ, АФ называется фильтром-прототипом, или просто прототипом.

Основными ограничениями для синтеза ЦФ по прототипу являются:

– сохранение существенных АЧХ прототипа в АЧХ ЦФ, что означает необходимость отображения мнимой оси p-плоскости на единичную окружность z-области;

– устойчивый прототип должен быть преобразован в устойчивый ЦФ, что означает необходимость отображения полюсов устойчивого прототипа из левой p-полуплоскости внутрь единичного круга z-области.

БИХ - фильтр рассчитывается на основе аналоговых ФНЧ - прототипов по следующей методике:

1. На основании требований к ЦФ формулируются требования к соответствующему аналоговому фильтру-прототипу.

2. Для аналогового фильтра-прототипа определяется нормированный аналоговый ФНЧ-прототип с частотой среза (границей полосы пропускания) и рассчитывается его нормированные полюсы и нули, исходя из квадрата его АЧХ.

3. С помощью формул преобразования частот производится преобразование нормированного ФНЧ-прототипа в аналоговый прототип, соответствующий исходному ЦФ.

4. Производится пересчет нулей и полюсов из аналоговой области в цифровую.

Идеальная (прямоугольная) форма АЧХ не может быть физически реализована, поэтому в теории аналоговых фильтров разработан ряд методов аппроксимации прямоугольных АЧХ. Аппроксимацией фильтра будем называть передаточную функцию, у которой АЧХ приближается к одной из идеальных характеристик. Такая передаточная функция должна удовлетворять следующим условиям:

1) она должна быть рациональной функцией от p с вещественными коэффициентами;

2) ее полюсы (нули знаменателя) должны лежать в левой полуплоскости p-плоскости;

3) степень полинома числителя должна быть равной или меньше степени полинома знаменателя.

Согласно методике расчета БИХ-фильтра, требования, предъявляемые у ЦФ, перекладываются на соответствующий аналоговый ФНЧ-прототип. Рассчитав ФНЧ, можно несложными преобразованиями изменить его частоту среза, превратить его в ФВЧ, полосовой, либо режекторный фильтр с заданными параметрами. Вследствие чего расчет аналогового фильтра начинается с расчета так называемого фильтра-прототипа, представляющего собой ФНЧ с частотой среза, равной 1 рад/с.

Рассчитаем АФ-прототип с передаточной функцией по требованиям, заданным к ЦФ, используя функции 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 0 0 0],

Передаточная функция АФ-прототипа

.

Частотные характеристики АФ-прототипа

Нули и полюсы АФ-прототипа

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

[r, p, k] = residue(b, a);

b, a – коэффициенты передаточной функции;

r вектор-столбец вычетов;

p вектор-столбец полюсов;

k - вектор-строка целой части дробно-рациональной функции.

Определим ИХ каждого звена, соответствующего :

.

Проведем их дискретизацию, заменой , получив при этом импульсные характеристики звеньев ЦФ.

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

Для вычисления коэффициентов полиномов воспользуемся функцией

p = ploy(A), где p – коэффициенты характеристического полином, A – его корни.

az = poly(exp(p*1/Fd));

az – коэффициенты знаменателя передаточной функции;

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]

Построим диаграмму полюсов и нулей

Увеличенная диаграмма нулей в области большего числа нулей

Все полюсы фильтра лежат внутри окружности единичного радиуса z-плоскости, следовательно, фильтр устойчив.

Составим разностное уравнение ЦФ и решим его методом итерационной процедуры для управляющего сигнала типа единичный скачок.

Решение методом итерационной процедуры в 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 – вектор управляющего сигнала типа единичный скачок.

Решение разностного уравнения для управляющего сигнала типа единичный скачок

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

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

Подав на вход фильтра с тестовый сигнал с единичной амплитудой и с частотой 1500 Гц, видим, что входной сигнал подавлен.

Подав на вход сигнал с частотой 7500 Гц, видим, что входной сигнал не подавлен, т. к. находится в полосе пропускания.

Подав на вход фильтра сигнал с частотой 14000 Гц, видим, что входной сигнал подавлен.

Из графиков видно, что в полосе задерживания сигнал практически полностью подавляется.

Составим каноническую структурную схемы фильтра. Для этого запишем уравнения:

где b0= ; b1= … b5= ;

a0=1 ; a1= … a6= .

Составим параллельную структурную схему фильтра. Для этого передаточную функцию представим в виде суммы простых дробей с помощью matlab-функции residue.

[r, p, k] = residue(bz, az);

bz, az – коэффициенты передаточной функции;

r вектор-столбец вычетов;

p вектор-столбец полюсов;

k - вектор-строка целой части дробно-рациональной функции.

; ; .

Проектирование и анализ фильтров в графической среде FDATool

Рабочее окно в FDATool (БИХ-фильтр)

АЧХ и ФЧХ БИХ-фильтра

Импульсная характеристика

Переходная характеристика

Все характеристики, построенные в FDATool для БИХ-фильтра совпадают с рассчитанными ранее.

Рабочее окно в FDATool (КИХ-фильтр)

АЧХ и ФЧХ КИХ-фильтра

Импульсная характеристика

Переходная характеристика

Диаграмма нулей

Все характеристики, построенные в FDATool для КИХ-фильтра совпадают с рассчитанными ранее.

Листинг 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;

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');

%--- Окно Кайзера

%-- определим 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]);

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--');

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)], [],'r--');

legend('АЧХ АФ-прототипа','АЧХ идеального фильтра'); title('АЧХ');

xlabel('f, Hz'); ylabel('Magnituda');

subplot(2,1,2); plot(w/(2*pi),unwrap(angle(h))*180/pi); grid;

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)];

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');

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');

Заключение

В данной курсовой работе был произведен расчет КИХ - и БИХ - фильтров. Расчет фильтров осуществлялся как стандартными средствами Matlab, так и с использованием пакета Signal Processing, обеспечивающий широкие возможности по созданию программ обработки сигналов для современных научных и технических приложений. Результаты расчетов идентичны, следовательно, синтез фильтров выполнен правильно.

Список использованной литературы

1. , Лазарева цифровых фильтров: Руководство к курсовой работе. – Чебоксары: Изд-во Чуваш. ун-т., 20с.

2.  Сергиенко обработка сигналов. - СПб.: Питер, 20с.

3.  , Антонов алгоритмы цифровой обработки сигналов: Метод. Указания к лабораторным работам. – Чебоксары: Изд-во Чуваш. ун-т., 2003, 60с.