цос / Лаб_раб_ №4

ЛАБОРАТОРНАЯ РАБОТА № 4

МОДЕЛИРОВАНИЕ ОБРАБОТКИ СИГНАЛОВ ЛИНЕЙНОЙ ДИСКРЕТНОЙ СИСТЕМОЙ

Цель работы: ознакомление с основными характеристиками линейных дискретных систем, изучение реакций системы на воздействия во временной и частотной областях

Необходимые теоретические сведения

Линейная система может быть представлена в нескольких эквивалентных формах, каждая из которых полностью описывает его:

• в форме рациональной передаточной функции (tf - представление); если система является непрерывной (аналоговой), то она описывается непрерывной передаточной функцией:

, (1)

где s- оператор Лапласа, s=jω.

Дискретная система может быть представлена в z-области дискретной передаточной функцией вида:

. (2)

В обоих случаях для задания звена достаточно задать два вектора – вектор b коэффициентов числителя num и а - знаменателя передаточной функции den.

• в виде разложения передаточной функции на простые дроби; в случае

простых корней такое разложение имеет вид (для дискретной передаточной функции):

(3)

где K - коэффициент усиления;

q_i - вещественный или комплексный нуль (корень числителя);

p_k - вещественный или комплексный полюс (корень знаменателя);

(m-1),(n-1) - количество нулей и полюсов соответственно.

• путем задания векторов z нулей передаточной функции, p - ее полюсов и k - коэффициента передачи звена (zp - представление):

(4)

• в каскадной форме (sos-представление).

• в пространстве состояний (ss- представление)

Пакет SignalTool предоставляет пользователю ряд процедур, позволяющих преобразовать систему из одной формы в другую. Если, например, известны полюса p₁=-0.2+0.5i и p₂=-0.2-0.5i, а также нули z₁=0.1i и z₂=-0.1i.

p=[-0.2+0.5i -0.2-0.5i]’;

z=[0.1i -0.1i]’;

k=5;

Здесь k - коэффициент, стоящий перед всей дробью передаточной функции (его можно выбрать произвольно, если нет иных соображений). Символ ‘ (апостроф), поставленный после закрывающей квадратной скобки, означает транспонирование. При этом вектор-строка преобразуется в вектор-столбец. (принято полюсы и нули задавать вектор-столбцами, а коэффициенты передаточной функции вектор-строками).

Преобразование полюсов и нулей в коэффициенты передаточной функции производится при помощи функции zp2tf (z,p,k):

[b,a]= zp2tf (z,p,k);

plot (p,’x’)

hold on

plot (z,’o’)

hold off

axis equal

axis ([-1 1 -1 1])

Таким образом можно получить диаграмму нулей и полюсов. Последняя команда позволяет задать граничные значения для действительной (первые два числа) и мнимой (следующие два числа) осей.

Для изображения карты нулей и полюсов можно использовать функцию zplane. Например, для заданных p и q имеем

Получить карту нулей и полюсов можно, указав коэффициенты числителя b (num) и знаменателя a (den) в порядке убывания степеней, начиная с коэффициента при нулевой степени. Первый элемент вектора den всегда равен 1. Например, так

Если требуется вывести значения каких-либо переменных рабочей области MatLab на монитор, то нужно просто после значка >> ввести имя

переменной, завершив ввод нажатием клавиши <Enter>. Например, если указать в командной строке z, то будет выведен вектор нулей. Вслед за оператором p будет выведен вектор полюсов. Если не ставить точку с запятой после операторов преобразования, то выходные параметры будут выведены на монитор без дополнительных команд.

Например, задана передаточная функция вида:

то записав следующие команды

имеем

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

Нахождение коэффициентов передаточной функции по коэффициентам

разложения ее на простые дроби (3) осуществляется использованием функций residue и residuez. Первая применяется для непрерывной передаточной функции вида (1), вторая - для дискретной передаточной функции (2). При обращении вида:

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

[ b, a] = residuez( r, p, k)

вычисляются коэффициенты числителя и знаменателя передаточной функции по заданным векторам ее разложения - вычетов r, полюсов p и коэффициентам целой части k.

С помощью тех же процедур осуществляется разложение заданной передаточной функции на простые дроби.

При этом обращение к ним должно быть таковым:

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

[ r, p, k] = residuez( b, a).

Для представления передаточной функции в виде суммы дробей (5)

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

Обратная процедура выполняется так:

Математическое моделирование обработки сигналов линейной системой (ЛС) включает:

• расчет реакции ЛС по передаточной функции;

• расчет характеристик ЛС во временной и частотной областях;

• анализ воздействия и реакции во временной и частотной областях.

Во временной области основными характеристиками ЛС является импульсная характеристика h (реакция на единичный импульс) и переходная характеристики g (реакция на единичный скачок).

В частотной области ЛС основной характеристикой является частотная характеристика:

,

а также ее модуль (АЧХ) и аргумент (ФЧХ).

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

Для аналоговой системы передаточная функция задается в форме tf -полиномов числителя и знаменателя: sys=tf(b,a). Например, так можно задать колебательное звено

kzv1 = tf([1 4], [1 2 100])

Transfer function:

s + 4

---------------

s^2 + 2 s + 100

Для дискретной системы задается значение Тs: sysd=tf(b,a,Ts).

Например,

kzv2 = tf([1 4], [1 2 100],0.01)

Transfer function:

z + 4

---------------

z^2 + 2 z + 100

Sampling time: 0.01

Систему, заданную как непрерывная, можно перевести в дискретную форму,воспользовавшись процедурой c2d:

sysd = c2d(sys, Ts, method)

Здесь sys - исходная непрерывная заданная модель, sysd - получаемая в результате работы процедуры дискретный аналог исходной системы, Ts - задаваемое значение шага дискретизации, method - параметр, определяющий метод дискретизации. Последний параметр может быть указан

'zoh' - соответствует применению экстраполятора нулевого порядка:

внутри интервала дискретизации сигналы аппроксимируются постоянной величиной, равной значению сигнала в начале интервала дискретизации;

'foh'- соответствует применению экстраполятора первого порядка: внутри интервала дискретизации сигналы аппроксимируются отрезками прямых,проходящих через концы кривой сигнала в интервале дискретизации;

'tustin' - билинейная аппроксимация Тастина внутри интервала дискретизации;

'prevarp' - та же аппроксимация Тастина с заданной частотой предыскривления;

'matched' - метод согласования нуля и полюса.

Процедура d2c осуществляет обратную операцию - переводит дискретную систему в непрерывную, например:

k1 = d2c(KZVd1)

Transfer function:

s + 4

---------------

s^2 + 2 s + 100

k2 = d2c(KZVd4,'tustin')

Transfer function:

s + 4

---------------

s^2 + 2 s + 100

Как можно убедиться, указанные операции являются взаимно-обратными. Процедура d2d позволяет переопределить дискретную систему, либо меняя шаг дискретизации

sys1 = d2d(sys,Ts),

либо вводя групповые задержки Nd (целое, в количестве шагов дискретизации)

sys1 = d2d(sys,[],Nd).

Приведем примеры. Вначале изменим шаг дискретизации на Ts=0.1 для системы KZVd1:

kd1=d2d(KZVd1,0.1)

Transfer function:

0.09352 z - 0.06018

-----------------------

z^2 - 0.9854 z + 0.8187

Sampling time: 0.1 ,

а затем введем задержку по входу, равную 3 Ts. Получим:

kd2=d2d(k d1,[],3)

Transfer function:

0.09352 z - 0.06018

-----------------------------

z^5 - 0.9854 z^4 + 0.8187 z^3

Sampling time: 0.1

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

impulse - нахождение отклика системы на единичное импульсное входное воздействие;

step - нахождение реакции системы на единичный скачок входного воздействия;

lsim - определение реакции системы на входное воздействие произвольной формы, задаваемое в виде вектора его значений во времени.

bode - строит графики АЧХ и ФЧХ (диаграмму Боде) указанной системы;

nyquist - строит в комплексной плоскости график Амплитудно-Фазовой

Характеристики (АФХ) системы в полярных координатах(годограф);

Например,

tf([0 -2 3 4 5],[5 8 0 0 -4])

step(tf([0 -2 3 4 5],[5 8 0 0 -4]))

impulse(tf([0 -2 3 4 5],[5 8 0 0 -4]))

nyquist(tf([0 -2 3 4 5],[5 8 0 0 -4]))

bode(tf([0 -2 3 4 5],[5 8 0 0 -4]))

Для применения процедуры lsim необходимо предварительно задать вектор 't' значений времени, в которых будут заданы значения входного воздействия, а затем и задать соответствующий вектор 'u' значений входной величины в указанные моменты времени.

t = 0:0.01:40; u = sin(t);

ssys=tf([0 -2 3 4 5],[5 8 0 0 -4])

lsim(ssys,u,t);grid

Для построения графиков АЧХ и ФЧХ для аналоговой системы нужно ввести векторы коэффициентов b и a в, а затем вызвать функцию расчёта комплексной частотной характеристики и построения графиков freqs (b,a).

Например, для b0=2, b1= -1, a0=1, a1=3, a2=2.5.

b=[2 –1];

a=[1 3 2.5];

freqs (b,a)

Будут построены графики АЧХ и ФЧХ (АЧХ в логарифмическом масштабе, но без пересчёта в децибелы, ФЧХ в градусах). По умолчанию выбираются 200 частот, логарифмически равномерно распределённых в диапазоне от 0.1 до 10. Если нужно построить АЧХ в линейном масштабе, в ином диапазоне частот, ФЧХ в радианах и с устранением скачков на 2πk радиан, то следует ввести следующие операторы (пусть, для примера, границы частотного диапазона от f1=100 Гц до f2=1000 Гц, число расчётных точек N =451):

f=100:(1000-100)/450:1000;

b=[1 0 7.971e6];

a=[1 7.427e2 1.501e6 5.536e8];

k=freqs (b,a,f*2*pi);

subplot (2,1,1)

plot (f,abs(k)/max(abs(k))), grid

subplot (2,1,2)

plot (f,unwrap(angle(k))), grid

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

impz - для расчета или графического отображения импульсной характеристики дискретной системы;

freqz – для расчета комплексного коэффициента передачи или построения графиков АЧХ и ФЧХ дискретной системы;

dbode – дискретный аналог диаграмм Боде;

grpdelay – для расчета группового времени задержки.

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

где N – количество точек, в которых рассчитывается импульсная характеристика;

Fs- частота дискретизации в Гц;

h- вектор-столбец отсчетов импульсной характеристики;

nT – вектор-столбец значений дискретного времени.

Можно вычислить импульсную характеристику, задав количество точек N:

Тогда

Для вычисления частотной характеристики H(e^jωt) по коэффициентам b (num) и a (den) используется функция freqz, формат может быть представлен в нескольких видах

где N – количество точек, в которых рассчитывается частотная характеристика;

w – вектор частот в радианах в секунду;

f- вектор частот в Гц;

Первый элемент den всегда равен 1.

Функция freqz является дискретным аналогом функции freqs. Например,

Те же значения можно получить, используя программу

Расчет АЧХ и ФЧХ дискретной системы можно выполнить с помощью функций abs и angle.

C помощью функции dbode , которая имеет формат, также можно вычислить АЧХ и ФЧХ:

где w– вектор значений частот ω от 0 до π/T

T- период дискретизации;

MAG - вектор значений АЧХ;

PHASE- вектор значений ФЧХ.

Например,

Групповое время задержки (ГВЗ) – это производная от ФЧХ. Расчет выполняется с помощью функции grpdelay:

Для предыдущего примера имеем:

Представленная дискретная передаточная функция описывает в сжатой форме такое конечно-разностное уравнение :

где b –вектор коэффициентов воздействия

a -вектор коэффициентов реакции

Моделирование дискретных систем можно призводить на основе уравнения свертки с нулевыми начальными условиями. Формула свертки имеет вид:

где импульсная характеристика h(m) и входное воздействие (x(n) задаются в виде конечных последовательностей .

В MatLab свертка задается при помощи функции соnv, формат которой имеет вид: conv(x,h) или conv(h,x)

где х - вектор воздействия длиной

- вектор отсчетов импульсной характеристики длиной

Результатом вычислений есть реакция системы длиной

Например, система задана разностным уравнением

В данной программе заданы только значения вектора b, вектор а=[1].

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

Функция deconv выполняет операцию, обратную свертке. Поэтому, если известна реакция (вектор у) и воздействие (вектор х), но не известны векторы коэффициентов b и a , импульсную характеристику можно найти как:

,

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

Необходимо помнить, что вычисление возможно, если первые элементы х и у не равны нулю! В рассмотренном примере это невозможно.

А для предыдущего примера x(n)=cos(0.5n) имеем:

Выполним обратную процедуру:

Как видно, коэффициенты b и h совпадают.

ЗАДАНИЕ

1. Согласно варианту задать систему в различных формах представления передаточной функции.

2. Записать разностное уравнение.

3. Построить карту нулей и полюсов.

4. Определить временные характеристики .

5. Найти частотные характеристики (в линейном и логарифмическом масштабах).

Вар.№	Полином числителя	Полином знаменателя
1	0.1	(1.5s*+1)(3s+1)^2
2	1.0	(((1/5)s)^2+2*1/5*0.1s+1)(3s+1)
3	5.0	(((1/5)s)^2+2*1/5*0.1s+1)s
4	10.0	(2.5s+1)(5s+1)^2
5	20.0	(((1/50)s)^2+2*1/50*0.1s+1)(5s+1)
6	30.0	(((1/50)s)^2+2*1/50*0.1s+1)s
7	3.0	(0.1s+1) ^2 (0.3s+1)
8	5.0	((5s)^2+2*5*0.1s+1) (5s+1)
9	7.0	(((1/50)s)^2+2*1/50*0.1s+1)s^2
10	s+1	(1.5s+1)(3s+1)^2
11	(s+2)(s-1)	(((1/5)s)^2+2*1/5*0.1s+1)(3s+1)
12	(s+5)(s-2)	(((1/5)s)^2+2*1/5*0.1s+1)s
13	(s+1)(s-1)	(2.5s+1)(5s+1)^2
14	((1/30s)^2+2*1/30*0.1s+1)	(((1/50)s)^2+2*1/50*0.1s+1)(5s+1)
15	1.0	(((1/50)s)^2+2*1/50*0.1s+1)s
16	5.0	(0.1s+1) ^2 (0.3s+1)
17	((3s)^2+2*3*0.1s+1)	((5s)^2+2*5*0.1s+1) (5s+1)
18	(s+1)^2	(((1/50)s)^2+2*1/50*0.1s+1)s^2
19	(1.5s+1)(3s+1)	(((1/50)s)^2+2*1/50*0.1s+1)s^2
20	8.0	(((1/50)s)^2+2*1/50*0.4s+1)(5s+1)

* s – переменная Лапласа.