i-808190579

														0.361																	0.323
									2			0.398															2T				0.282
														1.279																		1.469
														0.506																		0.67
3. Процедура вычисления коэффициентов рекуррентным МНК
		1015																					Звено 1
				A		I		1		0
								0		1
						C0 I
						1r0					0
											0
						for		k 1 N
												U
							k							k 1
											Y1
														k 1
							k								Ck 1 k
																T
											1				k				Ck 1 k
							C C										k			T C
								k				k 1							k				k 1
							e	k	Y1							T			1r			k 1
												k			k
							1rk 1rk 1												k ek
						1r
																					0.567												0.567
				i 0 N											A	N														1T
																					0.717												0.717
	0.6																									0.8
Ai	0.4																								Ai	0.6
																										1
	0
																											0.4
b1	0.2																								a1
																										0.2
	0				20			40				60					80				100					0					20		40		60		80		100
										i																								i

91

Звено 2

A2

I identity(4)

C1 I

0

0

2 1

0

0

for k 2 N

Y1

k 1

Y1k 2

k

Y2k 1

Y2

k 2

k

Ck 1 k

T

1

k

Ck 1 k

C C

k

T C

k

k 1

k k 1

e

k

Y2

T

2

k 1

k

k

2 k 2 k 1 k ek

2

0.361

A2N

0.398

1.28

0.506

92

			0.323
	2T		0.282
			1.469
			0.67

	1						4
A2i 0	0.5						A2i 1 3
	0	20	40	60	80	100	2
b21							b22
0.5							1
	1						0	20	40	60	80	100
				i						i

	8						10
	6
A2i 2 4						A2i	0	20	40	60	80	100
							3
a21	2						10
						a22
	0	20	40	60	80	100	20
	2						30
				i						i

Рис . П.4.2. Изменение оценок параметров звена 2

93

Тестовый пример в системе Matlab

clc;clear;close all; W=10; N=100; Ts=1;

u = random('Normal',W,W,1,N); V = normrnd(0,0.1*W,N,1); a1=0.7165; b1=0.5669;

a2 = [1.469 -0.6703];

b2 = [0.323 0.2824];

y1 = zeros(1,N);

y2 = zeros(1,N); y = zeros(1,N); for k=3:N

y1(k)=b1*u(k-1)+a1*y1(k-1); y2(k)=b2(1)*y1(k-1)+b2(2)*y1(k-2)+a2(1)*y2(k-1)+a2(2)*y2(k-2); t(k)=k;

y(k)=y2(k)+V(k);

end

data1 =iddata(y1',u',Ts); data4 = iddata(y2',u',Ts); plot(data4)

ident

94

Opening System Identification Tool ....... done. >> arx221

Discrete-time IDPOLY model: A(q)y(t) = B(q)u(t) + e(t) A(q) = 1 - 1.469 q^-1 + 0.6703 q^-2

B(q) = 0.323 q^-1 + 0.2824 q^-2 Estimated using ARX on data set W2

Loss function 2.88015e-028 and FPE 3.11056e-028 Sampling interval: 1

96

		ЛАБОРАТОРНАЯ РАБОТА № 5
МЕТОДЫ СГЛАЖИВАНИЯ И ФИЛЬТРАЦИИ
Цель работы – применение методов сглаживания и фильтрации времен-
ных рядов, оценка их эффективности.
	КРАТКОЕ ТЕОРЕТИЧЕСКОЕ ОПИСАНИЕ
Теоретические сведения приведены в разделах 4.6, 2.8, 2.9, 2.10 [1] и [2,
3].
Наблюдаемый временной ряд (рис. 5.1)можно представить в виде
		y(k) x(k) V (k) ,			k , , ...		(5.1)
где y(k) – наблюдаемое (измеряемое) значение временного ряда на k-ом ша-
ге; x(k) – регулярная (полезная) составляющая наблюдаемого временного
ряда на k-ом шаге; V (k) – аддитивная помеха на k-ом шаге наблюдения с ну-
левым математическим ожиданием m и конечной дисперсией .
			v			v
y(k)					y(N+m)
			x(k)
				y(k)
T0		T0
					mT0
0	1	2	k	N	Тп-интервал	k
	T-интервал наблюдения				упреждения
		Рис. 5.1. Наблюдаемый временной ряд

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

Функции сглаживания и фильтрации Mathcad.

Рассмотрим некоторые функции, используемые в Mathcad и приложениях Signal Processing Extension Pack, Data Analysis Extension Pack.

97

Метод экспоненциального сглаживания, expsmooth (y,α), реализует следующую процедуру

S(k) α y(k) (1 α)S(k 1) ,

k 1,2,...N , (5.2)

где S(k) , S(k ) – сглаженное значение временного ряда на k-ом и (k-1)-ом шаге наблюдения; y(k) – наблюденное (измеренное) значение временного

ряда на k-ом шаге наблюдения; – постоянная сглаживания; N – объем выборки.

Метод скользящей медианы

S:=medsmooth(y, N0)

возвращает вектор S сглаженных значений, размерности N . Здесь y – вектор наблюденных (измеренных) значений, размерности N , N – ширина

окна, по которому происходит сглаживание, оно должно быть нечетным числом, меньшим N. Сглаживание происходит в три этапа: сглаживается исходный вектор y процедурой medsmooth, далее вычисляются остатки, как разность исходного вектора y и сглаженного, вектор остатков сглаживается процедурой medsmooth и итоговый результат процедуры представляет сумму двух сглаженных векторов.

Адаптивное сглаживание обеспечивает подстройку параметров алгоритма сглаживания на каждом такте сглаживания, в зависимости от постулируемой модели тренда x(k) .

Процедуры приведены в разделе 4.6 [1].

Симметричная линейная процедура сглаживания методом наименьших квадратов по правилу N ближайших соседей

S:=supsmooth(x, y),

где x – вектор значений аргумента в возрастающем порядке размерности N , y – вектор измеренных значений временного ряда, размерности N . Параметр N процедура формирует адаптивно (различную ширину полосы

сглаживания для различных частей данных).

Фильтрация на основе дискретного преобразования Фурье.

Основное назначение одномерного дискретного прямого преобразования Фурье – выделить частоты регулярных составляющих сигнала, зашумленного помехами.

Если длина N последовательности наблюдаемых данных y(k) является

степенью числа 2, т.е. N 2n , то применяется алгоритм быстрого преобразования Фурье (БПФ), имеющий максимальную производительность.

Чтобы исключить влияние постоянной составляющей случайного процесса на картину спектра (сдвиг ноль-линии), необходимо его централизовать

Y0(k)=y(k)-yср.

98

Для фильтрации помех и выделения низкочастотного полезного сигнала x(k) вводят пороговую фильтрацию на полученные коэффициенты с j . В

спектре остаются частоты,

удовлетворяющие условию

m

c j

m , иначе

c j .

Здесь c

k

f

, частота на k-ом шаге такта, m

, m

– нижняя и верхняя

k

c

N

граница порога.

Применяя к отфильтрованному спектру обратное БПФ, получим оценку детерминированной основы x(k) .

Фильтрация на основе вейвлет-преобразования.

Волновое (вейвлетное) преобразование для ряда функций эффективнее БПФ. Если БПФ раскладывает исходную функцию в гармонический ряд, то при вейвлет-преобразовании используют специальные базисные функции, локализованные в некоторой ограниченной области аргумента.

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

В Mathcad использована функция Даубечи (Добеши). Для фильтрации помех также используют пороговую фильтрацию.

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

		N
SS		Sk xk ,	(5.3)
	N k

где Sk , xk – сглаженное и искомое значение временного ряда на k-ом шаге. Более подробно о функциях сглаживания и фильтрации в руководствах поль-

зователей Mathcad по ссылке Help→Functions→Curve fitting and smooth на английском языке, Help→Электронные книги→Data Analysis Extension Pack: Exponential Smoothing, а также Help→Электронные книги→Signal Processing→here→Index of Functions или на сайте [4] на русском языке.

Функции сглаживания и фильтрации в Matlab.

Функции сглаживания представлены в разделе Curve Fitting Toolbox. Метод скользящего среднего

S=smooth (y, N0).

Размер окна по умолчанию N0 5 .

Для обращения к другим методам сглаживания

S=smooth (y,’method’).

Сглаживание можно реализовать в диалоговом окне Data на вкладке Smooth. Окно Data открывается нажатием на кнопку Data в основном окне приложения сftool.

99

Для перехода в графический интерфейс в командном окне набираем команду, либо в m-файле

>>сftool.

Враскрывающемся списке Original data set выбираются данные, созданные в рабочей среде Matlab.

Встроке ввода Smoothed data set вводится имя, которое будет присвоено сглаженным данным.

Враскрывающемся списке Method выбирается один из способов сглаживания.

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

Для фильтра Савицкого-Голея требуется задание степени полинома для локальной регрессии, степень задается в строке ввода Degree.

Сглаженные данные можно визуализировать (кнопка View), переименовать (кнопка Rename), удалить (кнопка Delete) и сохранить в глобальной переменной рабочей среды Matlab.

С описанием остальных функций сглаживания, экспортом и импортом исходных данных, можно ознакомиться по ссылке Curve Fitting Tool-

box→Interactive Curve Fitting→(Data Fitting).

Фильтрация и сглаживание в приложении Signal Processing Toolbox (Digital Filter Implementation).

Функция

S=filter (b, a, y),

реализует дискретный фильтр

S(k) b0 y(k) b1 y(k 1) ... bm y(k m) a1S(k 1) an S(k n),

где n – порядок фильтра, b b0 b1...bm , a 1 a1...an – векторы коэффициен-

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

W (z) b0 b1z 1 ... bm z m . 1 a1z 1 ... an z n

Так для метода экспоненциального сглаживания, согласно разностного

уравнения (5.2)

b α , a 1 1 α , или a 1 α 1 .

Фильтр скользящей медианы

S=medfilt1(y,N0).

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

100