i-808190579

More Display Op-

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

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

c=fft(y-mean(y)).

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

Одномерное дискретное вейвлет-преобразование

[sa,sd]=dwt(y,’wname’),

реализует процедуру одноуровнего дискретного разложения данных, заданных вектором y, по базису вейвлет-функций wname, а также возвращает векторы, содержащие соответственно коэффициенты аппроксимации sa и детализации sd данного разложения. В качестве базисной функции используем функцию Добеши 4-го порядка.

wname=’db4’.

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

Обратное преобразование (восстановление данных) осуществляется функцией

X=idwt(sа, sd, ‘wname’);

возвращающей вектор коэффициентов аппроксимации Х для векторов sa и sd.

Удаление коэффициентов детализации sd x=idwt(sa,[ ],’wname’).

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

wavemenu.

В открывшемся окне Wavelet Toolbox Main Menu выбираем тип вейвлет преобразования Wavelet 1-D, в появившемся окне активизируем кнопку

File→Import from Workspace→Import→Signal→Y→OK. В окне Wavelet 1-D в

окне Loaded Signal отображается временной ряд y(t) . В правой половине ос-

новного окна выберем из меню Wavelet→db и рядом порядок 4, период квантования T0 1.

Далее активизируем кнопку Analyze. В левой половине появятся графики коэффициентов sa и sd по уровням.

Переключателем Display mode режимы вывода можно изменять, детали вывода можно устанавливать в отдельном подокне (кнопка

tions).

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

Окна Statistics и Histograms предназначены для анализа и графического

101

вывода статистических характеристик сигнала и всех коэффициентов его разложения. В окнах Compress и De-noise устанавливаются режимы компрессии (сжатия) сигналов и очистки сигналов от шумов.

ЗАДАНИЕ И ПОРЯДОК ВЫПОЛНЕНИЯ РАБОТЫ

Согласовать с преподавателем среду моделирования Mathcad или Mat-

lab.

Непрерывная модель сигнала и параметры приведены в табл. 5.1. Помеха распределена по нормальному закону с параметрами mv ,

V W .

Таблица 5.1

Параметры модели

Номер		Параметры модели				Детерминированная основа
варианта	а		b
1	5		-		-					x(t) a
2	1		1		-			x(t) a b t
3	3		- /4		1	x(t) a sin t b
4	4		0.04		-		x(t) a aebt
							x(t)			b		a
							a									t
							0				6·a
5	1		5		-				T :		a
												N
						X (t) : if a t 6 a,b,0
							t : ,T ... a
6	2		0.1	v	0.1 W	x(t)						a
					T 1
										a exp( bt)
					0
7	7		-		-					x(t) a
8	2		0.5		-			x(t) a b t
9	1		0		0.5	x(t) a sin t b
10	10		0.02				x(t) a aebt

1. Выберите в соответствии с вариантом W модель x(t) , и согласно

тестового примера разработайте программу расчетов.

Для T . с. осуществите переход к дискретной модели и рассчитайте дискретный временной ряд x(kT0 ) объемом N n , для n . Наблюдаемый временной ряд

102

yk xk Vk , k ,...,N .

Постройте график x(k) и y(k) . Рассчитайте статистики центрированного временного ряда Ryy (r) , S yy ( ) , постройте их графики и сформулируйте суждение о структуре модели x(k) .

2. Реализуйте сглаживание временного ряда y(k) рассмотренными

методами и постройте график x(k) и сглаженных значений ~ . yk

Рассчитайте эффективность методов по критерию (5.3).

3. Реализуйте дискретное преобразование Фурье и постройте график спектральной плотности с j . Исходя из графика выберите уровни m , m по-

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

~

ний y(k) реализуйте обратное преобразование Фурье. Постройте график x(k) и восставленных значений временного ряда yf (k) . Меняя порог m , до-

бейтесь наилучшего качества фильтрации по критерию (5.3).

4.Пункт 3 повторите для вейвлетного преобразования.

5.Повторите задание для Vk .

6.Для среды Matlab дополнительно получить результаты в графическом режиме GUI.

Выводы по работе должны включать:

1)обоснование выбора интервала квантования T ;

2)анализ Ryy (r) и S yy ( ) ;

3)анализ эффективности методов сглаживания и выбор наиболее эффективного. Соответствие модели и метода сглаживания;

4)влияние параметров сглаживания для каждого метода на эффективность сглаживания;

5)влияние уровня помех V на эффективность сглаживания.

КОНТРОЛЬНЫЕ ВОПРОСЫ

1.Необходимость сглаживания временных рядов.

2.Что представляет собой детерминированная основа наблюдаемого временного ряда?

3.Какие методы можно применять для сглаживания временных ря-

дов?

4.Охарактеризуйте процедуру выбора модели детерминированной

основы.

5.Рекомендации по выбору начальных условий процедур сглажива-

ния.

103

6.Рекомендации по выбору параметров сглаживания.

7.От чего зависит точность прогнозирования в задачах прогнозирования временных рядов?

8.Можно ли применить МРА к сглаженным значениям временного

ряда?

ЛИТЕРАТУРА

1.Масальский, Г.Б. Математические основы кибернетики. Часть 1. Основы идентификации систем управления: учебное пособие [Электронный ресурс] / сост. Г.Б. Масальский. – Красноярск: Сиб. федер. ун-т, 2014. – 1 диск.

2.Макаров Е. Инженерные расчеты в Mathcad 15: Учебный курс. –

СПб.: Питер, 2011. – 400 с.: ил.

3.Дьяконов В.П. Matlab. Полный самоучитель. – М.: ДМК Пресс,

2012, - 768 с.: ил.

4.Образовательный математический сайт Exponenta.ru. Режим досту-

па: http://www.exponenta.ru/

104

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

1. Модель временного ряда

W 2

n 7

N 2n 128

v W

V rnorm(N 0 v )

a 2

b 4

1

T0 0.1

X(t) a sin( t b)

f0

1 10

0 2 f0

62.832

T0

k 0 N 1

x

a sin( k T0 b)

y

k

x

V

ycp mean(y) 0.161

k

k

k

y0 y

k

ycp

k

tk k

- вектор значений аргумента для процедуры fft

10

5

xk

0

yk

5

100

50

100

150

k

N r 1

y0k y0k r

r 0

N

Ryyr

1

Ryy0 5.918

2

N r

k 0

1

RY correl (y0 y0)

RY0

757.497

Ryr N r

RYr

Ry0 5.918

Py cfft (Ryy)

Py0 1.65

Sy cfft (y0)

Sy0 0

length(Py) 65

SY fft (y0)

SY0 0

length(Sy) 128

length(SY) 65

105

	6						15
	4					P yr
							10
Ryyr 2						Syr
Ryr	0					SYr	5
	2
	4	0	20	40	60	80	0
								0	2	4	6	8	10
				r						r
2. Сглаживание
				0.3		N0 5
S expsmooth(y )						S1 medsmooth(y N0)
S0 Var(x S) 0.705					SS1 Var(x S1) 1.308
	4
xk	2
Sk	0
S1k	2
	4
	60				50			100					150
						k

	4
	2
xk	0
medsmooth(y 3)k 2
	4
	60	50	100	150
			k

106

3. Адаптивное сглаживание

yy1

a00

y0

yy

a0

y0

a1

0

0

for

k 1 N 1

for

k 1 N 1

yyk ak 1

yy1

k

a0

k 1

a1

k 1

(1 ) yy

a

k

y

k

k

y

2

yy1

k

a0

y

(1

)

y

k

k

k

yy

k

2 yy1

a1

k

a1

k 1

k

y

k

y1 supsmooth (t y)

yy1

ss Var(x yy)

0.76

ss1 Var(x yy1) 2.222

ss2 Var(x y1)

0.374

10

xk

5

yyk

0

y1k

yy1

k

5

100

50

100

150

k

Сглаживание с пос ледующей с плайн-интерполяцией

z cspline(t y)

X(t) interp(z t y1 t)

4

xk

2

X(t)

0

20

50

100

150

k t

107

4. Дискретное преобразование Фурье

C fft (y)

j 0 N

m2 max

C

19.605

m1 5

20

N

Cj

15

i 0

2

10

ydi if m1

m2 Ci0

m1

Ci

5

yf ifft (yd)

00		20	40	60
			j
					Ss Var(x yf) 0.068
		3
		2
xk		1
yfk		0
	1
	2
	30			50	100	150
					k
5. Вейвлетное преобразование

	w wave (y)	m3 6	m4	20
15		j 0 N 1
wk 10		swj if m3 wj m4 wj 0
m3		yw iwave(sw)
5
		Sw Var(x yw)		0.232

00	50	100	150
		k
	4
	2
xk
ywk	0
	2
	40		50	100	150
				k

108

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

clc;clear;scrsz = get(0,'ScreenSize');
n=7;
N=2^n;	% объѐм выборки
W=2;	% номер варианта
%
k = 1:N;	% счѐтчик
x = W*sin(0.05*k+pi/4);	% детерминированная основа
cv =W;	% СКО помехи
V = random('Normal', 0 , cv, 1, N);	% реализация закона распределения по-
мехи
y = x+V;	% Реализация временного ряда

%% Отображение графика временного ряда

hFig = figure('MenuBar','none','NumberTitle','off','Name','Характеристики временного ряда:',...

'Units', 'pixels', 'OuterPosition', [0 0 scrsz(3) scrsz(4)]); subplot(2,1,1);

plot(k,y,k,x); grid on;

legend('Временной ряд','Детерминированная основа');

%%Построение графиков корреляционной и спектральной функций

Ryy = xcov(y,'unbiased'); C = fft(y);

r = 1:(N/2);

%Отображение графиков корреляционной функции subplot(2,2,3);

plot(r,Ryy(r)); grid on;

legend('Автокорреляционная функция Ryy(r)'); %Отображение графиков спектральной функции subplot(2,2,4);

plot(r,abs(C(r))); grid on;

legend('Спектр-Фурье fft(y) ');

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

N0=5;	% Ширина окна сглаживания
S = medfilt1(y,N0);
%% Экспоненциальный фильтр-------------------------------------------------
alfa=0.3;		% постоянная сглаживания (вес)
epha = filter(alfa, [1 alfa-1], y);		% реализация экспоненциального
фильтра
%Метод скользящего среднего--------------------------------------------------
-----
oa = smooth(y,N0);
% Графики	скользящего среднего, экспоненциального и скользящей медианы
hFig1 = figure('MenuBar','none','NumberTitle','off','Name','Графики
фильтров:',...
'Units', 'pixels', 'OuterPosition', [0 0 scrsz(3) scrsz(4)]);
subplot(2,1,1);
plot(k,S,k,epha,k,oa,k,x);

grid on;

legend('Скользящей медианы','Экспоненциальный','Скользящего среднего','Детерминированная основа');

%% Дискретное преобразование Фурье----------------------------------------

m1=70;	%Уровень отсечения высоких частот
for i=1:N
if	m1<=abs(C(i));
	109

yf(i)=C(i); else

yf(i)=0;

end

end yf=ifft(yf);

%%Одномерное дискретное вейвлет-преобразование-----------------------------

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

[sa,sd]= dwt(y,'db4'); yd = idwt(sa,[],'db4');

%%Отображение графика Вейвлет преобразования и преобразования Фурье

subplot (212);

plot (k,yf,k,yd,k,x); grid on;

legend('Фурье','Вейвлет','Детерминированная основа');

hFig2 = figure('MenuBar','none','NumberTitle','off','Name','Графики фильтров:',...

'Units', 'pixels', 'OuterPosition', [50 100 500 200]); %% Отображение Эффективности

axis('off');

hFig2=text(0.0,1,' Квадрат невязки:','FontSize',16); hFig2=text(0,0.8,sprintf('Скользящей медианы = %g',var(x-S)),'FontSize',16); hFig2=text(0,0.6,sprintf('Экспоненциального = %g',var(x- epha)),'FontSize',16);

hFig2=text(0,0.4,sprintf('Скользящего среднего = %g',var(x- oa')),'FontSize',16);

hFig2=text(0,0.2,sprintf('Фурье = %g',var(x-yf)),'FontSize',16); hFig2=text(0,0,sprintf('Вейвлет = %g',var(x-yd)),'FontSize',16);

110