Лекция 10. Статистически анализ

ЛЕКЦИЯ 10. СТАТИСТИЧЕСКИ АНАЛИЗ

Очевидно ВКФ се явява усредненото по време значение на произведението на първата функция умножена на втората функция, която е отместена във времето със време на задръжка .

И така взаимната спектрална плътност (ВСП) на два стационарни процеса се задава от формулата:

(2)
Сега можем да преминем към дискретния вариант на горните формули когато процесите са зададени в ограничен интервал от време Т и в n на брой дискретни точки разположени равноотстоящи на интервал T_s. Тогава формула (1) се трансформира в следния вид:

(3)
или в по-упростен вид:

(4)
Дискретизираната формула (2) има следния вид:

(7)
от което следва, че ВСП на два процеса, за всяка стойност на честотата, е равна на произведението на комплексния спектър на втория процес умножен на комплексно спрегнатото Фурие-изображение на първия процес на същата честота.

С помощта на филтър може да се формира случаен процес със зададена корелационна функция. Отначало генерираме нормално разпределен /Гаусов случаен процес/, след което пропускаме сигнала през филтър, чийто динамически свойства определят корелационната функция на резултатния случаен сигнал.

1. 
   Ще формираме два случайни Гаусови процеса x1(t) и x2(t), които се
   

Първият случаен процес формираме с командите както следва:

Ts=0.01; t=0:Ts:20; x1=randn(1,length(t));

plot(t,x1);grid

Вторият процес е следният:

Ts=0.001; t=0:Ts:20; x2=randn(1,length(t));

plot(t,x2);grid

2. 
   Ще създадем дискретен филтър от втори порядък със собствена честота
   

Ts=0.01; t=0:Ts:20; x1=randn(1,length(t));omega0=2*pi; zeta =0.05; A=1; omegas=omega0*Ts;

a=[1+2*zeta*omegas+omegas^2, -2*(1+zeta*omegas), 1];

b=[A*2*zeta*omegas^2];

y1=filter(b,a,x1);

plot(t,y1); grid

За вторият процес само времевата стъпка е по-малка на порядък т.е. имаме:

Ts=0.001; t=0:Ts:20; x2=randn(1,length(t));omega0=2*pi; zeta =0.05; A=1; omegas=omega0*Ts;

a=[1+2*zeta*omegas+omegas^2, -2*(1+zeta*omegas), 1];

b=[A*2*zeta*omegas^2];

y2=filter(b,a,x2);

plot(t,y2); grid

Вижда се от получените графики, че и в двата случая на изхода на формиращия филтър се образуват случайни процеси с преобладаваща честота 1 Нz.

Фурие преобразование на случаен процес

Ще направим Фурие преобразование на първия случаен процес дефиниран в предишния параграф във вид на бял Гаусов шум със времева стъпка 0.01 и продължителност 100 с.

Ts=0.01; T=100; t=0:Ts:100; x1=randn(1,length(t));omega0=2*pi; zeta =0.05; A=1; omegas=omega0*Ts;

a=[1+2*zeta*omegas+omegas^2, -2*(1+zeta*omegas), 1];

b=[A*2*zeta*omegas^2];

y1=filter(b,a,x1);

plot(t,y1); grid

Ще изчислим Фурие-изображението (ФИ) на формирания случаен шум

като отчетем вече разгледаната процедура за отстраняване на преходните процеси и ще построим графика на модулите на ФИ и на СПМ.

%формираме честотните масиви

df=1/T; Fmax=1/Ts; f=-Fmax/2:df:Fmax/2;dovg=length(f);

%изчисляваме коригираните масиви на Фурие-изображенията

Fu1=fft(x1)/dovg; Fu2=fft(y1)/dovg;

Fu1p=fftshift(Fu1); Fu2p = fftshift(Fu2);

%изчисляваме модулите на Фурие-изображенията

A1=abs(Fu1p); A2 = abs(Fu2p);

%изчисляваме спектралната плътност

S1=Fu1p.*conj(Fu1p)*dovg; S2=Fu2p.*conj(Fu2p)*dovg;

%графики на бял шум

subplot(2,1,1); stem(f,A1); grid; subplot(2,1,2);stem(f,S1),grid

%графики на филтрирания нормален /Гаусов/ шум

c1=fix(dovg/2)-200; c2=fix(dovg/2)+200;

subplot(2,1,1); stem(f(c1:c2), A2(c1:c2)); grid;

subplot(2,1,2); stem(f(c1:c2), S2(c1:c2)); grid

В разгледания пример използвахме връзката между спектрална плътност и Фурие-изображение за да определим спектралната плътност на случаен процес. Обаче в MatLab (Signal Processing Toolbox) има процедура, която позволява директно да се определи СП на сигналите. Това може да стане с командата:

[S, f]=psd(x, nfft, Fmax);

Името на процедурата е абревиатура от psd =power spectral density, x е вектора, който задава поредицата от стойности на изследвания процес; nfft е броя на елементите на вектора x, които се обработват с процедурата fft; Fmax=1/Ts е честотата на дискретизация на сигнала; S е вектора който задава стойностите на спектралната плътност; f е вектора на честотите, за които са определени стойностите на спектъра. В общия случай дължината на последните два вектора е равна на nfft/2.

Ще демонстрираме процедурата psd() за намиране спектъра на плътноста на същия случаен процес както в предишния пример, с доминираща честота 1 Нz.

%формираме случаен процес със нормално разпределение

Ts=0.01; T=100; t=0:Ts:100; x1=randn(1,length(t));omega0=2*pi; zeta =0.05; A=1; omegas=omega0*Ts;

a=[1+2*zeta*omegas+omegas^2, -2*(1+zeta*omegas), 1];

b=[A*2*zeta*omegas^2];

y1=filter(b,a,x1);

%формираме честотния вектор

df=1/T; Fmax=1/Ts; f=-Fmax/2:df:Fmax/2;dovg=length(f);

%прилагаме процедурата psd()

[S,f]=psd(y1,dovg, Fmax);

%извеждаме резултат графично

stem(f(1:200), S(1:200)), grid

Ако приложим същата процедура без да присвоим върнатия резултат на изходен вектор автоматично ще се изведе графика на СП, като значенията на спектъра се извеждат в логаритмичен мащаб в децибели. Пробвайте следната команда на командната линия на МatLab:

psd(y1, dovg, Fmax)

Корелационни функции

Процедурата xcorr() позволява да се изчисляват както взаимната корелационна функция (ВКФ ) на два процеса така и автокорелационната функция (АКФ) на един процес. Ако x и у са векторите с дължина n на два процеса то командата c=xcorr(x,y); позволява да се изчисли вектора c на ВКФ с дължина 2n-1.

Нека като пример изчислим АКФ на разглеждания вече случаен процес y1:

R=xcorr(y1); tau=-10+Ts : Ts : 10; lt=length(tau);

s1r=round(length(R)/2)-lt/2;

s2r=round(length(R)/2)+lt/2-1;