Преобразование Лапласа и z-преобразование — конспект лекции

Преобразование Лапласа

Преобразование Лапласа — это интегральное преобразование, связывающее функцию F(s) комплексного переменного (изображения) с функцией f(x) вещественного переменного (оригиналом). Используется для решения дифференциальных и интегральных уравнений, а также для анализа динамических систем. Нас сейчас интересует последнее.

Прямое преобразование Лапласа для вещественного f(t) выглядит так:

(1)   \begin{equation*} F(s) = \int\limits_{0}^{\infty} f(t) e^{-st} \,dt, \end{equation*}

где:

  • f(t) – оригинал;
  • F(s) – изображение функции f(t);
  • s = \sigma + j \omega;
  • \sigma — некое число;
  • \omega — комплексная частота, рад./с.

Т.к. e^{-st} не имеет размерности, s должен иметь размерность “1/время”, или размерность частоты. Поэтому переменную Лапласа s называют комплексной частотой. А ещё e^{-st} представляет собой общую форму решения линейных дифуров.

По сути, преобразование Лапласа можно рассматривать как непрерывную функцию, значение которой при некотором s представляет собой корреляцию функции f(t) и затухающей синусоиды e^{-st}, частота которой равна \omega, а коэффициент затухания \sigma. Давайте запишем e^{-st} следующим образом:

(2)   \begin{equation*} e^{-st} = e^{-\sigma t - j \omega t} = \frac{e^{-j \omega t}}{e^{\sigma t}} \end{equation*}

e^{-j \omega t} — вектор, модуль которого равен единице, вращающийся вокруг начала координат в комплексной плоскости с частотой \omega. Ещё его называют фазором. e^{\sigma t} — комплексное число, равное единице при t=0, значение которого увеличивается с увеличением t. Таким образом, e^{-st} на комплексной плоскости представляет собой спираль:

Спираль, описывающая множитель в преобразовании Лапласа
Функция e^{-st} при \omega = 2\pi, \sigma = 1

Ниже показана действительная часть e^{-j \omega t} при разных значениях \omega, и \sigma:

Действительная часть множителя преобразования Лапласа
Действительная часть e^{-st} при разных значениях \omega, и \sigma

Получается, что в процессе расчёта F(s) мы делаем корреляцию нашего сигнала с сигналом, имеющим различную сходимость и частоту колебаний. Отсюда можно сделать вывод о поведении анализируемой нами системы. Одной из основных характеристик системы является устойчивость. Простыми словами: система является устойчивой, если при ограниченном входном сигнале, на выходе также получается ограниченный сигнал, который после устранения входного воздействия самостоятельно возвращается к некоторому установившемуся значению.

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

(3)   \begin{equation*} H(s) = \frac{Y(s)}{X(s)} \end{equation*}

Давйте разберём пример. Возьмём передаточную функцию некоторой системы, запишем её изображение в виде:

(4)   \begin{equation*} H(s) = \frac{s^3+2s^2-s}{s^2+3s+2} \end{equation*}

И построим график |H(s)|:

b = [1 2 -1 0];
a = [1 3 2];

omega = linspace(-2.5, 2.5);
sigma = linspace(-2.5, 0.5);

[sigmagrid, omegagrid] = meshgrid(sigma, omega);
sgrid = sigmagrid + 1i*omegagrid;
H = polyval(b, sgrid)./polyval(a, sgrid);

mesh(sigma, omega, abs(H));
xlabel('σ');
ylabel('jω');
zlabel('|H(s)|');

Результат выполнения данного скрипта:

Функция |H(s)|
Функция |H(s)|

Мы наблюдаем 3D-поверхность, устремляющуюся в бесконечность в двух точках: при j\omega=0, \sigma = -2 и j\omega=0, \sigma = -1. Это полюсы системы. Для того, чтобы найти полюсы системы, нужно прировнять её знаменатель к нулю и найти корни получившегося уравнения. Если проделать аналогичную процедуру с числителем, т.е. найти точки, где H(s)=0, мы получим нули системы.

Если взять и посмотреть сечение поверхности |H(s)| при \sigma = 0, получим ни что иное, как преобразование Фурье данной системы, или её АЧХ. Это легко проверить, подставив в уравнение (1) s=j\omega:

(5)   \begin{equation*} F(s) = \int\limits_{0}^{\infty} f(t) e^{-j\omega t} \,dt \end{equation*}

Правда, знакомая формула? Да это и есть преобразование Фурье непрерывного сигнала.

Давайте подробнее рассмотрим полюсы системы. Для этого построим их в двумерной плоскости, где по оси x будет \sigma, по оси yj\omega. В Matlab это делается очень просто. Добавим в листинг 1.1 пару строчек:

Hs = tf(b, a);         % передаточная функция
iopzmap(Hs), grid on;  % график нулей и полюсов

Получим следующий результат:

Нули и полюсы системы H(s)
Нули и полюсы системы H(s)

Нули отмечены кружками, полюсы — крестиками.

Если все полюсы системы расположены слева от оси j\omega, система является устойчивой. Если хоть один полюс расположен справа от этой оси — система является неустойчивой. Если все полюсы расположены на оси j\omega, система находится на границе устойчивости (условно устойчива), к таким системам относят, например, генераторы.

Z-преобразование

Преобразование Лапласа применяют для непрерывных систем, а для анализа интересующих нас дискретных систем используют z-преобразование. Для дискретного сигнала x[n] z-преобразование X[z] выглядит так:

(6)   \begin{equation*} X[z] = \sum\limits_{-\infty}^{+\infty} x[n] \cdot z^{-n}, \end{equation*}

где z — комплексное число. Запишем его как z = re^{j\omega}, где r — модуль, а \omega — аргумент комплексной переменной.

Получается, при r=1 выражение (6) принимает вид:

(7)   \begin{equation*} X[z]|_{z=e^{j\omega}} = \sum\limits_{-\infty}^{+\infty} x[n] \cdot e^{-j\omega n} \end{equation*}

Выражение (7) представляет из себя ДПФ от сигнала x[n]. В общем случае, X[z] представляет собой поверхность, если взять её сечение цилиндром |z|=1, то на поверхности данного цилиндра будет отображаться АЧХ нашего дискретного сигнала.

Анализ устойчивости здесь производится аналогично анализу с помощью преобразованию Лапласа. За исключением того, что мы смотрим расположение полюсов не относительно оси j\omega, а относительно единичной окружности |z|=1, центр которой находится в начале координат. Если все полюсы находятся внутри единичной окружности, система является устойчивой. Если хоть один полюс находится снаружи — система неустойчива. Полюсы, расположенные на единичной окружности говорят об условной устойчивости.

Рассмотрим пример в пакете Matlab. Возьмём систему с передаточной характеристикой:

(8)   \begin{equation*} G(s) = \frac{0.1s^2+10s+5}{s^3+8s^2+4s+2} \end{equation*}

Это будет прямая связь. И вторую систему:

(9)   \begin{equation*} C(s) = \frac{2s+1}{2s+3} \end{equation*}

Это будет обратная связь. Далее дискретизируем их с помощью функции c2d (в результате чего получаем Gd[z] и Cd[z] соответственно) и соединяем согласно схеме, показанной на рисунке ниже (в нашем случае K=1):

Структурная схема анализируемой системы
Структурная схема анализируемой системы

Затем строим график нулей и полюсов получившейся системы на комплексной плоскости:

clear

G = tf([0.1 10 5],[1 8 4 2]);
Gd = c2d(G,0.1);
C = tf([2 1],[2 3]);
Cd = c2d(C,0.1);
sys = feedback(Gd,Cd);

figure
rlocus(sys), grid on;

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

Корневой годограф системы из листинга
Корневой годограф системы из листинга “Анализ устойчивости дискретной системы”

Но, помимо этого, на графике видны кривые разных цветов, которые показывают траекторию движения нулей и полюсов при разном коэффициенте усиления K. Получается, что мы с вами построили корневой годограф. С его помощью можно отследить, при каких значениях K система будет устойчивой, а при каких — нет, а также величину пререгулирования при ступенчатом воздействии на сигнал (становится видно, если нажать левой кнопкой мыши на графике). В нашем случае граница устойчивости системы — K=21. При большем коэффициенте усиления система становится неустойчивой.

Построим графики реакции системы из предыдущего листинга на ступенчатое воздействие (функция Хевисайда, или ступенька — это сигнал, который при t<0 равен нулю, при t\geq 0 равен единице). Дополним код:

figure
subplot(2,2,1)
sys1 = feedback(Gd,Cd);
step(sys1), grid on
title('K=1')

subplot(2,2,2)
sys2 = feedback(Gd,Cd*10);
step(sys2), grid on
title('K=10')

subplot(2,2,3)
sys3 = feedback(Gd,Cd*20);
step(sys3), grid on
title('K=20')

subplot(2,2,4)
sys4 = feedback(Gd,Cd*25);
step(sys4), grid on
title('K=25')

Ниже показан результат выполнения данного кода:

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

Из рисунка видно, что с увеличением K появляются колебания на фронте переходной характеристики. При K=25 (больше 21) система не возвращается в состояние равновесия (она неустойчива).

Вернёмся к БИХ-фильтрам. Помните, в его структурной схеме были прямоугольники с надписью “Задержка”? Давайте посмотрим, как она выглядит в z-области. Рассмотрим на примере задержки на 1 такт:

(10)   \begin{equation*} y[n] = x[n-1] \end{equation*}

Запишем z-преобразование для уравнения (10):

(11)   \begin{equation*} Y[z] = \sum\limits_{-\infty}^{+\infty} x[n-1]z^{-n} \end{equation*}

Пусть k=n-1, тогда:

(12)   \begin{equation*} Y[z] = \sum\limits_{-\infty}^{+\infty} x[k] z^{-(k+1)} = \sum\limits_{-\infty}^{+\infty} x[k] z^{-k}z^{-1} = z^{-1}\sum\limits_{-\infty}^{+\infty} x[k] z^{-k} \end{equation*}

или:

(13)   \begin{equation*} Y[z] = X[z] \cdot z^{-1} \end{equation*}

Получается, чтобы сделать задержку на 1 такт, достаточно домножить сигнал на z^{-1}, поэтому во многих структурных схемах фильтров вместо надписи “Задержка” можно увидеть просто z^{-1}. Задержка на n тактов выглядит как z^{-n}.

Выражение для БИХ-фильтра в z-области будет выглядеть следующим образом:

(14)   \begin{equation*} Y[z]=X[z]\sum\limits_{k=0}^{N} b[k] z^{-k}+Y[z]\sum\limits_{k=1}^{M} a[k] z^{-k} \end{equation*}

Далее, чтобы составить передаточную характеристику, всё, что относится к выходному сигналу, отнесём в числитель, а входного — в знаменатель. Получим:

(15)   \begin{equation*} H[z] = \frac{Y[z]}{X[z]} = \frac{\sum\limits_{k=0}^{N} b[k] z^{-k}}{1-\sum\limits_{k=1}^{M} a[k] z^{-k}} \end{equation*}

Другим не менее важным параметром фильтра является его частотная характеристика. Чтобы её найти, достаточно в его передаточную функцию (15) подставить z=e^{j\omega}:

(16)   \begin{equation*} H[\omega] = H[z]|_{z=e^{j\omega}} = \frac{\sum\limits_{k=0}^{N} b[k] e^{-jk\omega}}{1-\sum\limits_{k=1}^{M} a[k] e^{-jk\omega}} \end{equation*}

Пока это всё.


Скачать конспект в pdf: Laplace and z-transform Lecture – V.V. Leonidov.pdf