Спектрална плътност на мощноста
Ще дефинираме понятието спектрална плътност на мощноста (СПМ) или съкратено просто спектрална плътност (СП). Това понятие се дефинира като Фурие преобразование на корелационната функция R12() и се прилага обикновено за два едновременно протичащи стационарни процеса x1(t) и x2(t). Взаимната корелационна функция (ВКФ)на два такива процеса се определя от формулата:
(1)
Очевидно ВКФ се явява усредненото по време значение на произведението на първата функция умножена на втората функция, която е отместена във времето със време на задръжка .
И така взаимната спектрална плътност (ВСП) на два стационарни процеса се задава от формулата:
(2)
Сега можем да преминем към дискретния вариант на горните формули когато процесите са зададени в ограничен интервал от време Т и в n на брой дискретни точки разположени равноотстоящи на интервал Ts. Тогава формула (1) се трансформира в следния вид:
(3)
или в по-упростен вид:
(4)
Дискретизираната формула (2) има следния вид:
(5)
Ако сега заместим (4) в (5) и променим реда на сумирането можем да получим следната връзка между ВСП и резултатите /означени с у/ от процедурата fft:
(6)
където чертата отгоре означава комплексно спрегната стойност. Формула (6) може да се запише и по следния начин:
(7)
от което следва, че ВСП на два процеса, за всяка стойност на честотата, е равна на произведението на комплексния спектър на втория процес умножен на комплексно спрегнатото Фурие-изображение на първия процес на същата честота.
ФОРМИРАНЕ НА СЛУЧАЕН ПРОЦЕС
С помощта на филтър може да се формира случаен процес със зададена корелационна функция. Отначало генерираме нормално разпределен /Гаусов случаен процес/, след което пропускаме сигнала през филтър, чийто динамически свойства определят корелационната функция на резултатния случаен сигнал.
-
Ще формираме два случайни Гаусови процеса x1(t) и x2(t), които се
отличават по времевата стъпка Ts, т.е по гъстотата на случайните пулсации.
Първият случаен процес формираме с командите както следва:
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
-
Ще създадем дискретен филтър от втори порядък със собствена честота
1Hz т.е. w0=2 и коефициент на затихване =0.05. За първия процес имаме:
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;
plot(tau, R(s1r:sr2)); grid
Сподели с приятели: |