i-808190579

Кxx:=correl(x,x),

где Кxx – вектор, размерности вектора х.

В приложении Matlab Signal Processing Toolbox расчет ведется по фор-

муле (2.5).

Обращение к функции

Кxx=xcorr(x,’unbiased’).

Оценка автокорреляционной функции центрированного стационарного временного ряда равна

	^						N r
	R xx (r)						x(k) x(k r) , r			, , ,...,m . (2.6)
				N		r
							k
	В приложении Mathcad:
	covar (x,x)r= xk mx xk r mx ;
					k
			Rxx:=covar (x,x)
	и в Matlab
	Rxx=xcov(x,’unbiased’).
	Оценка нормированной корреляционной функции равна
			^							(2.7)
			(r) R xx (r) , r , , ,...,m .
	^
		xx
			^
					R xx ( )

В приложении Mathcad

ρxx:=lcorr(x,x),

где ρхх – вектор размерности вектора х. В Matlab соответственно

rхх=xcov(x,’coeff’).

Если имеются реализации центрированных стационарных временных рядов u(k) и x(k) , то оценки их взаимных корреляционных функций рассчи-

тываются по формуле (2.5) и (2.6) соответствующей заменой одной переменной x u :

Rux=xcov(u,x,’unbiased’).

Оценка спектральной плотности случайного процесса. Для оценки спектральной плотности стационарного процесса применяются непараметрические и параметрические методы (см. разд. 2.10 [1]).

В непараметрических (nonparametric) методах используются наблюдения временного ряда x(k) . В параметрических (parametric) методах исходят

из того, что наблюдаемый случайный процесс представляет собой результат воздействия дискретного «белого шума» на динамическое звено n-го порядка. Спектральный анализ сводится в данном случае к решению оптимизационной задачи, то есть поиску таких параметров фильтра, при которых выход фильтра y(k) наиболее близок к наблюдаемому сигналу x(k) . Для этого хо-

рошо подходят авторегрессионные методы. Эффективность оценок опреде-

21

Signal Processing

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

Оценка спектральной плотности в приложении

Mathcad.

Sxx(r):=pspectrum(x,d,α,[w]) – плотность спектра мощности, Sux(r):=cspectrum(u,x,N,r,[w]) – плотность кросс-спектра,

где x, u – вектор реализаций временных рядов x(k) и u(k); d – целое число в диапазоне 1<d<length(x), определяет число перекрывающихся сегментов временного ряда x(k) ; α – постоянная сглаживания 0<α<1; обычно α=0.5; w –

индекс, определяющий окно сглаживания первичных оценок (периодограммы), для окна Хэннинга w=4, при w=0 данные не сглаживаются.

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

Оценка спектральной плотности мощности в приложении Signal Processing Toolbox Matlab.

Вычисление периодограммы непосредственно по наблюдениям случайного процесса x(k) (непараметрический метод)

[Pxx,w]=periodogram(x,window,w),

где Pxx – вектор оценки спектральной плотности мощности (PSD) размерности N / 2 1;

w – вектор нормированной угловой частоты от 0…π;

window – вектор коэффициентов окна сглаживания (по умолчанию прямоугольное окно).

Для уменьшения пульсаций (сглаживания) периодограммы можно воспользоваться непараметрическим методом Уэлча (Welch)

[Pxx,w]=pwelch (x,window,noverlap,w);

где noverlap – доля перекрытия соседних сегментов, по умолчанию 50%.

В методе Уэлча вектор наблюдений x(k) делится на перекрывающиеся

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

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

Непараметрический метод оценки плотности кросс-спектра.

Pxy=cpsd(x,y,window,noverlap).

Обозначения аналогичны предыдущим функциям.

Параметрический метод оценки на основе авторегрессионного метода Берга (Burg)

[Pxx,f]=pburg(x,n,nfft,f0),

22

где n – порядок авторегрессионной модели, nfft – число точек N для БПФ, f0 – частота квантования.

Для апериодических процессов типа «цветной шум» n=1. Для случайных процессов с гармонической составляющей он должен быть в два раза больше числа предполагаемых синусоидальных колебаний.

Число точек для расчета спектра по умолчанию nfft=256. Выходные возвращаемые векторы по умолчанию имеют длину nfft и соответствуют

полному диапазону частот 0…f0 для комплексного вектора x(k) и от 0...f20

для вещественного вектора x(k) .

Интерактивный режим спектрального анализа в Matlab [7].

В командном окне набираем >>sptool. Появившиеся три окна программы позволяют:

загрузить входные данные;

просмотреть их график;

осуществить их спектральный анализ предполагаемыми методами;

рассчитать фильтр;

провести фильтрацию сигнала рассчитанным фильтром.

Загрузка данных. Для загрузки сигнала в главном окне программы sptool активизируем File→Import.

Впоявившемся окне Sourse активизируем рабочую область Matlab (From Workspace). В окне Workspace Contents появятся переменные (векто-

ры) рабочей области Matlab.

Втретьем окне в поле Import AS выбираем вариант Signal.

Выберите в списке окна Workspace Contents переменную, содержащую отсчеты загружаемого сигнала и щелкните на кнопке - →, расположенной рядом с полем ввода Data. Можно и вручную ввести в это поле идентификатор переменной.

Вполе ввода Sampling Frequency по умолчанию находится значение

Fs=1. Ее можно отредактировать вручную, или импортировать, ввести идентификатор переменной или выбрав его в списке и щелкнув на нижней кнопке

-→.

Вполе ввода Name можно отредактировать имя, под которым данный сигнал (вектор введенных данных) будет фигурировать в программе SPTool (удобно оставить исходный идентификатор). Сделав это, щелкните по кнопке ok. Импортированный сигнал появится в списке Signals основного окна. Процедуру повторить для всех вводимых из рабочей области данных.

Просмотр графика. Выберем имя сигнала в списке Signals и щелкнем на кнопке View, расположенной под эти окном внизу. Появится окно Signal Browser. В окне Signals выбираем имя файла для просмотра.

23

Спектральный анализ. Выберем сигнал в списке Signals и активизируем кнопку Create, расположенной под списком (окном) Spectra. В появившемся окне Spectrum Viewer в левой части, выберем метод спектрального анализа. Данные для настройки параметров можно взять из проведенного ранее расчета, представленного в строке Inherit from, представленного в списке Spectra основного окна программы. После настройки параметров активизируем кнопку Apply. Будет рассчитана оценка спектра сигнала в соответствии с выбранным методом и выведен соответствующий график.

Расчет фильтра. Для расчета дискретного фильтра щелкните на кнопке New Design, расположенной под списком (окном) Filters в главном окне программы SPTool. Можно изменить параметры уже рассчитанного фильтра, выбрав его в списке Filters и щелкнув кнопку Edit Design. В появившемся окне Filter Designer необходимо задать все параметры для синтеза фильтра, руководствуясь специальной литературой, например [7].

Фильтрация сигнала. Для этого необходимо выбрать сигнал и фильтр соответственно в списках Signals и Filters основного окна программы и далее щелкнуть на кнопку Apply, расположенной под списком Filters. В появившемся окне Apply Filter в поле Output Signal задаем имя для выходного сигнала. В раскрывающемся списке Algorithm можно выбрать функцию Matlab для осуществления фильтрации – filter и другие. Выполнив эти действия, щелкните на кнопке ok. Будет рассчитан выходной сигнал, который появится под указанным именем в списке Signals основного окна программы.

Просмотреть график выходного сигнала и выполнить анализ его спектра можно описанным ранее способом. Для одновременного просмотра графика входного и отфильтрованного (выходного) сигнала необходимо их выбрать в списке Signals основного окна программы. Ввод каждого сигнала осуществляется с зажатой клавишей Ctrl и далее щелкаем на кнопке View.

Сохранение результатов. Сеанс работы можно сохранить с помощью команды Save Session или Save session As из меню File основного окна программы. Загрузить сохраненный сеанс можно командой Open Session того же меню File. Результаты расчетов можно экспортировать обратно в рабочее пространство Matlab командой Export из меню File. В появившемся окне Export from SPTool выбираются экспортируемые объекты и после щелчка на кнопке Export to Workspace, данные перемещаются в рабочую область

Matlab.

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

ритмы методов оценивания в Mathcad и Matlab. Для анализа различий следует обратиться к документации соответствующей функции.

24

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

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

lab.

1. Для среды Mathcad согласно тестового примера и вашего варианта получить, объяснить, проанализировать полученные результаты и сделать выводы. Параметры динамического звена K W , T W , где W – вариант,

N2n , n 9 , T0 0.25 Т .

2.Для среды Matlab согласно тестового примера и вашего варианта

получить, объяснить, проанализировать полученные результаты и сделать выводы. Из пакетного шагового режима перейти в интерактивный режим GUI и получить основные характеристики для процессов u(k) и x(k) . Сопос-

тавить и объяснить их с результатами пакетного режима. 3. Повторите исследования для T0 0.1T .

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

1) по пункту 1 задания объяснение случайного процесса, реализации случайного процесса, сечения случайного процесса для момента времени t и t , а также одномерной плотности вероятности; сопоставление полученных результатов с результатами примера 1.5 в [1], обоснование выбора интервала квантования T ; влияние дискретизации на полученные оцен-

ки; сравнительную оценку алгоритмов вычисления u , Su ;

2) по пункту 1, 2 задания анализ оценок корреляционных функций и спектральных плотностей временного ряда uk и временного ряда xk , закона

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

нительный анализ закона распределения выходного процесса нелинейного звена;

3) оценку влияния T на статистические характеристики временного ряда x(k) .

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

1.Определение случайного процесса и реализации случайного про-

цесса.

2.Определение случайной величины X t .

3.Какой процесс полностью характеризуется одномерной плотностью вероятностей?

4.Понятие «белого» и «цветного» шума.

25

5.Примеры стационарных и нестационарных случайных процессов.

6.Среднее значение по множеству и по времени, условия отличия и совпадения.

7.Для каких процессов корреляционная функция со временем стремиться к нулю?

8.Понятие времени корреляции.

9.Можно ли по корреляционной функции определить дисперсию случайного процесса?

10.Определение спектральной плотности.

11.Примеры спектральных плотностей.

12.Меняются ли статистические характеристики входного сигнала при прохождении через динамическое звено?

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

14.Изложите суть дискретизации непрерывного процесса.

15.Как выбрать интервал квантования T ?

16.Как выбрать граничную частоту c ?

17.Что такое дискретный временной ряд?

18.В чем преимущество и недостатки рекуррентных процедур вычисления выборочных средних и дисперсий?

19.Что характеризуют корреляционные функции?

20.Что характеризуют спектральные плотности?

ЛИТЕРАТУРА

1.Масальский, Г.Б. Математические основы кибернетики: учебное пособие с грифом СибРО УМО вузов РФ / Г.Б. Масальский. – Красноярск:

Сиб. федер. ун-т, 2012. – 176 с.

2.Теория автоматического управления: В 2-ч ч. Ч.2 / Под ред. А.В.

Нетушила. – М.: Высш. шк., 1972. – 432 с.

3.Отнес Р., Эноксон Л. Прикладной анализ временных рядов. Основ-

ные методы. – М.: Мир, 1982. – 428 с.

4.MATHCAD 6.0 PLUS. Финансовые, инженерные и научные расчеты

всреде Windows 95. Издание 2-е, стереотипное – М. Информационноиздательский дом «Филинъ», 1997 – 712 с.

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

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

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

2012. – 768 с.: ил.

7.Сергиенко А.Б. Цифровая обработка сигналов: Учебное пособие. –

СПб.: Питер, 2002. – 608 с.; ил.

26

				Тестовый пример в системе Mathcad
1. Моделирование случайного процесса U(k)
Ис ходные данные
W 3			вариант	mu 0	u 2	n 9	N 2n			N 512
T 4			K 3		T0 0.25 T		T0 1
Генератор с лучайных чис ел
u0 rnorm(N mu u )				u1 rnorm(N mu u )				u2 rnorm(N mu u )
Генератор шумов				v u gaussn (N)
k 0 N 1				uk u0k
	10						10
u1k	5						5
u2k	0					vk	0
u0k	5						5
	10	0	200	400	600		10	0	200	400	600
			k							k
				Рис . 2.1. Реализации с лучайных процес с ов
2. Статистики входного случайного процесса U(k)

Usr mean(u)

Usr

S2u Var(u)

S2u

- выборочная дис перс ия

Us0 u0

k 1 N 1

Su Stdev(u)

Su

СК О

Usk Usk 1

uk Usk 1

- рекуррентный алгоритм длявыборочной с редней

k

k 2 N 1

Us

u0

u1

Su

u

Us

2

u

Us 2

1

0

1

1

2

1

1

k 2

u

k

Us

k 1

2

Su

Su

- рекуррентный алгоритм длявыборочной

k 1

k

k 1

k

дис перс ии

27

	1											6
	0.4											4.8
Us r 0.2 0				200			400	600			S2u	3.6
Usk 0.8											Suk	2.4
1.4												1.2
	2											0	200		400	600
						k								k
Рис . 2.2.Изменение выборочного с реднего и дис перс ии при вычис лении рекуррентными
									алгоритмами
k 0 N 1							U0k uk Usr			- центрирование временного ряда
		N					- макс имальное чис ло шагов," макс имальное запаздывание"
m ceil
			4
r 0 m 1							- номер шага
			1			N r 1				- оценка автокорреляционной функции
Ruu(r)							U0k U0k r
		N								центрированного временного ряда
				r			k 0
Ruu(0) 3.999							Ru correl (U0 U0)				Ru 2.048 103
					1							0
	Ruc(r)						Rur		Ruc(0) 3.999
				N r
uu(r)			Ruu(r)				-оценка нормированной корреляционной функции
			Ruu(0)
	u lcorr(U0 U0)								u0 1
	4											1
	3											0.7
Ruu(r)	2											0.4
												ur
Ruc(r)	1											0.1
												0	30	60	90	120 150
	0	30			60		90	120	150			0.2
1												0.5
							r								r
		Рис . 2.3.Оценкикорреляционных функций временного ряда u(k)
											28

Rufr Ruu(r)

Su fft (Ruf)

Su0 0.178

SuU fft (U0)

SuR pspectrum (Ruf 1 0.1 0)

SuR0 0.038

SuU5 0.771 1.33i

Sup pspectrum (U0 8 0.5 0)

Sup0 3.717

r

length(Sup) 112

u

r

I 0 length(Sup)

m

J 0 length(SuR)

length(SuR) 127

1.5

15

Sur

1

2 SuU

10

r

SuRJ 0.5

SupI

5

00

0.1

0.2

0.3

0.4

00

0.2

0.4

0.6

0.8

1

J

ur lengthSuR(

)

ur

I

)

lengthSup(

Рис . 2.4. Оценка однос торонней с пектральной плотнос ти методом Б ПФ

3. Cтатистики дискретного случайного процесса X

k

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

	T0	b K (1	a)
a exp	T
		a 0.779

xk 1 a xk b uk

x0 K mu

x0 0

b 0.664

Xsr mean(x) Xsr 0.639 S2x Var(x)

S2x 4.577

- начальные ус ловия

-выборочное

среднее;

-выборочная дис перс ия

Оценка закона рас пределения
	p 1 3.2 log(N)		n ceil(p)		n 10
	int hh 0	h hh 1
	10
	5
xk	0
	5
	100	200	400	600
		k

hh histogram(n x)

150
100
h
50
0	5
10		0	5	10
		int

Рис . 2.5. Временной ряд x(k) и его гис трогамма

29

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

clc;

clear all;

%Исходные данные n=9;N=2^n;mu=0;su=2;T=4;K=3;T0=T*0.25;tm=N/4;

%Три реализации случайного процесса типа "белый шум":U0(t),U1(t),U2(t). for t=1:N

U(t)=normrnd(mu,su);%Для получения нормально распределенных случайных чи-

сел

U1(t)=normrnd(mu,su);%используем стандартную функцию matlab - normrnd. U2(t)=normrnd(mu,su);

tk(t) = t; %вектор времени (аргумента) для sptool

end

%Строим график t=1:N;

plot(t',[U' U1' U2'],'LineWidth',1); grid on

%Обозначаем оси xlabel('t');

title('Реализации случайного процесса (СП) белый шум(нормального распределения)'); %Добавляем легенду кривых

hleg1 = legend('U(t)','U1(t)','U2(t)'); pause

%Расчет и печать Kuu ,Ruu случайного процесса U(t)

[Kuu,ti]=xcov(U,'coeff'); [Ruu,ti]=xcov(U,'unbiased'); for j=1:tm

Ru(j)=Ruu(N-1+j); ru(j)=Kuu(N-1+j); tn(j)=j;

end plot(tn-1,Ru,'--',tn-1,ru); title('Оценка Ruu,ru'); hleg1 = legend('Ruu','ru'); grid on

pause

%Оценка и печать спектральной плотности мощности (СПМ) usr = mean(U);

SU2=var(U);

U0 = U-usr;

[Pu,w] = periodogram(U0,[],[],1); Pucp=mean(Pu);

[Pwu,ww] = pwelch(U0, [], [],[],1); Pwucp=mean(Pwu); plot(w,Pu,'--',ww,Pwu);

grid on

hleg1 = legend('Оценка СПМ функцией periodogram для U0' ,...

'и функцией pwelch для U0 '); hold off

pause

%Прохождение реализации случайного процесса через дискретный фильтр

a= exp(-T0/T);

b= (1-a)*K; x(1)=K*mu; for k=1:N-1

x(k+1)=a*x(k)+b*U(k);

end

30