Лазарев Ю.Ф. Mатематическое моделирование физических процессов и технических систем в MATLAB

Подождите немного. Документ загружается.

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

проксимации полиномов рис. 1.18:

x = -0.45:0.1:0.45;

y = [-1.1 0.2 0.1 0.8 0.5 0.2 0.4 0.1 0.1 -0.6];

stem(x,y); hold on

P1=polyfit(x,y,1)

P2=polyfit(x,y,2)

P3=polyfit(x,y,3)

P4=polyfit(x,y,4)

x1 = -0.5 : 0.01 : 0.5;

y1=polyval(P1,x1);

y2=polyval(P2,x1);

y3=polyval(P3,x1);

y4=polyval(P4,x1);

plot(x1,y1, '--',x1,y2,':',x1,y3, '.',x1,y4),

grid, set(gca, 'FontSize', 12),

title('Полиномиальная аппроксимация ');

xlabel('Аргумент'); ylabel('Функция');

legend('исходные','данные','Аппроксимация','линейная','квадратическая','кубическая','четв

ертой степени',0)

Рис. 1.18. Результаты применения функции polyfit

Функция

spline(X,Y,Xi) осуществляет интерполяцию кубическими сплайнами. При обращении

Yi = spline(X,Y,Xi)

она интерполирует значение вектора Y, заданного при значениях аргумента, представленных в векторе Х, и

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

Xi. В случае, если вектор Х не указан, по умолчанию принимается, что он имеет длину вектора Y и любой его

элемент равен номеру этого элемента.

В качестве

примера рассмотрим интерполяцию вектора

x = -0.45:0.1:0.45;

y = [-1.1 0.2 0.1 0.8 0.5 0.2 0.4 0.1 0.1 -0.6];

x1 = -0.46:0.01:0.46;

y2 = spline(x,y,x1);

stem(x,y); hold on

plot (x1,y2, '.'), grid

set(gca, 'FontSize',12),

title('Интерполяция процедурой SPLINE ');

xlabel('Аргумент'); ylabel('Функция')

Результат приведен на рис. 1.19.

Рис. 1.19. Интерполяция функцией spline

Одномерную табличную интерполяцию осуществляет процедура interp1. Обращение к ней в общем случае

имеет вид:

Yi = interp1(X,Y,Xi,’<метод>’),

и позволяет дополнительно указать метод интерполяции в четвертом входном аргументе:

'nearest'

ступенчатая интерполяция;

'linear'

линейная;

‘cubic'

кубическая;

‘spline'

кубическими сплайнами.

Если метод не указан, осуществляется по умолчанию линейная интерполяция. Например, (для того же вектора):

y = [-1.1 0.2 0.1 0.8 0.5 0.2 0.4 0.1 0.1 -0.6];

x1 = -0.46:0.01:0.46;

y1 = interp1(x,y,x1);

y4 = interp1(x,y,x1,'nearest');

y2 = interp1(x,y,x1,'cubic');

y3 = interp1(x,y,x1,'spline');

plot (x1,y1, '--',x1,y2,'.',x1,y3,x1,y4,':'), grid

set(gca, 'FontSize',12),

legend('линейная','кубическая','сплайновая','ступенчатая',0)

title('Интерполяция процедурой INTERP1 ');

xlabel('Аргумент'); ylabel('Функция')

Результаты приведены на рис.1.20.

Рис. 1.20. Интерполяция функцией interp1

1.4.5. Векторная фильтрация и спектральный анализ

В системе MatLAB есть несколько функций для проведения цифрового анализа данных наблюдений (измере-

ний).

Так, функция

y = filter(b,a,x) обеспечивает формирование вектора y по заданным векторам b, a, x в соот-

ветствии с соотношением:

y(k) = b(1)*x(k) + b(2)*x(k-1) + ... + b(nb+1)*x(k-nb) --a(2)*y(k-1) - a(3)*y(k-3) - ... - a(na+1)*y(k-na), (1.1)

где вектор b имеет такой состав

b = [ b(1), b(2),... , b(nb+1)],

a вектор а

a = [1, a(2), a(3),... , a(na+1)].

Соотношение (1) можно рассматривать как конечно-разностное уравнение фильтра с дискретной передаточной

функцией вида рациональной дроби, коэффициенты числителя которого образовывают вектор b, а знаменателя

вектор а, на вход которого подается сигнал x(t), а на выходе формируется сигнал y(t).

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

времени, соответствующие заданным значениям входного сигнала x(t) (вектор х).

Ниже приведен пример применения функции

filter

» x = 0:0.1:1;

» b = [1 2];

» a = [ 1 0.1 4];

» y = filter(b,a,x)

y =

0 0.1000 0.3900 0.2610 -0.5861 0.3146 3.9129

0.2503 -13.4768 2.8466 56.4225

Функции fft

(Fast Fourier Transformation) и

ifft (Inverse Fast Fourier Transformation) осуществляют преобра-

зования заданного вектора, соответствующие дискретному прямому и обратному преобразованиям Фурье.

Обращение к этим функциям вида:

y = fft ( x, n ); x = ifft ( y, n )

приводит к формированию вектора y в первом случае и х - во втором по формулам:

∑

−⋅−⋅⋅−

⋅=

nkmj

emxky

/)1()1(2

)()(

; (1.2)

∑

−⋅−⋅⋅

⋅=

nkmj

eky

/)1()1(2

)(

, (1.3)

где j - обозначение мнимой единицы; n - число элементов заданного вектора х (оно представляет также размер

выходного вектора y).

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

функции, являющейся суммой двух синусоид с частотами 5 и 12 Гц (рис. 1.21).

t =0:0.001:2;

x = sin(2*pi*5*t) + cos(2*pi*12*t);

figure, plot(t, x); grid,

set(gcf,'color','white'), set(gca,'FontName','Arial Cyr','FontSize',16),

title('Входной процесс '); xlabel('Время (с)'); ylabel('Х(t)')

Рис. 1.21. Бигармонический процесс

Найдем Фурье-изображение этого сигнала и выведем графическое представление модуля его Фурье-

изображения:

y = fft(x); a =abs(y);

plot(a); grid, set(gca,'FontName','Arial Cyr','FontSize',16),

title('Модуль Фурье - изображения ');

xlabel('Номер элемента вектора'); ylabel('abs(F(X(t))')

Результат отображен на рис. 1.22.

Рис. 1.22. Модуль Фурье-изображения

Теперь осуществим обратное преобразование с помощью функции

ifft

и результат также выведем в форме

графика:

z = ifft(y); plot(t, z);

grid, set(gca,'FontName','Arial Cyr','FontSize',16),

title('Обратное Фурье-преобразование ');

xlabel('Время (с)'); ylabel('Z(t)')

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

Рис. 1.23. Результат обратного преобразования Фурье

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

а) номер m отвечает моменту времени t

, в который измерен входной сигнал x(m); при этом ;

б) номер k - это индекс значения частоты f

, которому отвечает найденный элемент y(k) дискретного преобра-

зования Фурье;

в) чтобы перейти от индексов к временной или частотной области, необходимо знать значение h дискрета (ша-

га) времени, через который измерен входной сигнал x(t) и промежуток T времени, на протяжении которого он

измеряется; тогда шаг (дискрет) по частоте в

изображении Фурье определится соотношением:

Df = 1/T, (1.4)

а диапазон изменения частоты - формулой

F = 1/h; (1.5)

так, в анализируемом примере (h =0. 001, T = 2, n = 21):

Df = 0.5; F = 1000;

г) из (2) следует, что индексу k = 1 отвечает нулевое значение частоты (f

= 0); иначе говоря, первый элемент

вектора y(1) является значением Фурье-изображения при нулевой частоте, т. е. - просто суммой всех заданных

значений вектора x; отсюда получаем, что вектор y(k) содержит значение Фурье-изображения, начиная из час-

тоты f

= 0 (которой отвечает k = 1), до максимальной частоты fmax = F (которой отвечает k = n); таким обра-

зом, Фурье-изображение определяется функцией

fft

только для положительных частот в диапазоне от 0 до

F; это неудобно для построения графиков Фурье-изображения от частоты; более удобным и привычным пред-

ставляется переход к вектору Фурье-изображения, определенному в диапазоне частот от (-F/2) до F/2; часто-

та

получила название

2/FF

частоты Найквиста;

д) как известно, функция

является периодической по z с периодом 2 ; поэтому информация об Фурье-

изображении при отрицательных частотах расположена во второй половине вектора y(k).

⋅

Сформируем для анализируемого примера массив частот, исходя из вышесказанного:

f = 0 : 0.5 : 1000;

и выведем график с аргументом-частотой (рис. 1.24):

plot(f,a); grid, set(gca,'FontName','Arial Cyr','FontSize',16),

title('Модуль Фурье - изображения ');

xlabel('Частота (Гц)'); ylabel('abs(F(X(t))')

Рис. 1.24. Модуль Фурье-изображения в частотной области

Как следует из рассмотрения рис. 1.24, по нему непросто распознать те частоты (5 и 12 Гц), с которыми изменя-

ется входной сигнал. Это - следствие того обстоятельства, которое было отмечено в примечании г). Чтобы оп-

ределить частотный спектр входного сигнала, нужно сначала преобразовать полученный вектор y Фурье-

изображения с помощью процедуры

fftshift

Функция

fftshift (обращение к ней осуществляется таким образом: z = fftshift(y)) предназначена для

формирования нового вектора z из заданного вектора y путем перестановки второй половины вектора y в

первую половину вектора z . При этом вторая половина вектора z состоит из элементов первой половины векто-

ра y. Более точно эту операцию можно задать соотношениями:

z(1) = y(n/2+1); ... , z(k) = y(n/2+k);... , z(n/2) =

y(n); z(n/2+1) = y(1);...

..., z(n/2+k) = y(k);... z(n) = y(n/2).

Примечание. Операцию

fftshift удобно использовать для преобразования массива Фурье-

изображение с целью построения его графика в частотной области. Тем не менее,

этот массив не может быть использован для обратного преобразования Фурье.

Проиллюстрируем применение этой функции к предыдущему примеру (рис. 1.25):

f1 = -500 : 0.5 : 500; % Перестройка вектора частот

v = fftshift(y); % Перестройка вектора Фурье-изображения

a =abs(v); % Отыскание модуля

plot(f1(970:1030),a(970:1030)); grid, % Вывод графика

set(gca,'FontName','Arial Cyr','FontSize',16),

title('Модуль Фурье - изображения');

xlabel('Частота (Гц)'); ylabel('abs(F(X(t))')

Рис. 1.25. Преобразование Фурье-изображения функцией fftshift

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

тами 5 и 12 Гц.

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

избежание этого, нужно весь вектор y Фурье-изображения разделить на число его элементов

(n), чтобы полу-

чить вектор комплексного спектра сигнала:

N=length(y); a =abs(v)/N;

plot(f1(970:1030),a(970:1030)); grid

set(gca,'FontName','Arial Cyr','FontSize',14,'Color','white'),

title('Модуль комплексного спектра');

xlabel('Частота (Гц)');ylabel('abs(F(X(t)) / N')

Рис. 1.26. Комплексный амплитудный спектр сигнала

Результат приведен на рис. 1.26. Из его рассмотрения вытекает, что "амплитуды" всех составляющих гармоник

равны 0.5. При этом нужно принять во внимание, что "амплитуды" распределены между положительными и

отрицательными частотами поровну, поэтому они вдвое меньшее действительной амплитуды соответствующей

гармоники.

1.5. Построение простейших графиков

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

женерной и научной практике. MatLAB предоставляет для этого широкие возможности. Ознакомимся с важ-

нейшими из средств построения графиков, предусмотренными системой.

1.5.1. Процедура plot

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

режиме калькулятора.

Основной функцией, обеспечивающей построение графиков на экране дисплея, является функция

plot. Общая

форма обращения к ней такова:

plot(x1, y1, s1, x2, y2, s2,...)

Здесь x1, y1 - заданные векторы, элементами которых являются массивы значений аргумента (х1) и функции

(y1), отвечающие первой кривой графика; x2, y2 - массивы значений аргумента и функции второй кривой и т.д.

При этом предполагается, что значения аргумента откладываются вдоль горизонтальной оси графика, а значе-

ния функции - вдоль вертикальной оси. Переменные s1, s2,... являются символьными (их указание не

является

обязательным). Любая из них может содержать до трех специальных символов, определяющих соответственно:

а) тип линии, которая соединяет отдельные точки графика; б) тип точки графика; в) цвет линии. Если перемен-

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

в такой очередности: - синий

, зеленый, красный, голубой, фиолетовый, желтый, черный и белый - в зависимо-

сти от того, какая по очереди линия выводится на график. Например, обращение вида

plot(x1,y1,x2,y2,...) приведет к построению графика, в котором первая кривая будет линией из отрезков

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

Графики в MatLAB всегда выводятся в отдельное графическое окно, которое называют фигурой.

Приведем пример. Пусть нужно вывести график функции

y = 3sin(x +

/3)

на промежутке от -3

до +3 с шагом /100.

π π π

Сначала надо сформировать массив значений аргумента х:

x = -3*pi : pi/100 : 3*pi,

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

y = 3*sin(x+pi/3)

и, наконец, построить график зависимости y(х).

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

» x = -3*pi:pi/100:3*pi;

» y = 3*sin(x+pi/3);

» plot(x,y)

В результате на экране появится окно с графиком (см. рис. 1.27).

Рис. 1.27. График функции с указанием аргумента

Если вектор аргумента при обращении к функции

plot не указан явно, то система по умолчанию принимает в

качестве аргумента номера элементов вектора функции. Например, если ввести команду

» plot(y),

то результатом будет появление графика в виде, приведенном на рис. 1.28.

Графики, приведенные на рис. 1.27, 1.28, имеют следующие недостатки:

- на них не нанесена сетка из координатных линий, что затрудняет чтение графиков;

- нет общей информации о кривой графика (заголовка);

- неизвестно, какие величины отложены по осям графика.

Рис. 1.28. График функции без указания аргумента

Первый недостаток устраняется с помощью функции

grid

. Если к этой функции обратиться сразу после об-

ращения к функции

plot:

» x = -3*pi:pi/100:3*pi;

» y = 3*sin(x+pi/3);

» plot(x,y), grid,

то график будет снабжен координатной сеткой (рис. 1.29).

Ценной особенностью графиков, построенных в системе MatLAB, является то, что

сетка координат всегда

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

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

Рис. 1.29. Результат применения функции grid

Заголовок графика выводится с помощью процедуры

title

. Если после обращения к процедуре

plot

вызвать

title

таким образом:

title(‘<текст>’),

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

текст всегда должен помещаться в апострофы.