��Akademia G�rniczo-Hutnicza Katedra Robotyki i Mechatroniki Identyfikacja i analiza sygnaB
�w Laboratorium 10 Metody nieparametryczne w identyfikacji i analizie sygnaB
�w Przetwarzanie sygnaB
�w z zastosowaniem procedur statystycznych Signal Processing Toolbox zawiera zestaw wa|
niejszych funkcji do estymacji pewnych charakterystyk sygnaB
�w losowych. W szczeg�lno[
ci s
to narz
dzia do estymacji funkcji korelacji i kowariancji ci
gu danych oraz funkcje g
sto[
ci widmowej mocy sygnaB
�w dyskretnych. Korelacja i kowariancja Funkcje xcorr i xcov estymuj
wzajemn
korelacj
i kowariancj
ci
g�w proces�w losowych. Autokorelacja i autokowariancja s
traktowane jako specjalne przypadki korelacji i kowariancji wzajemnej. Korelacja wzajemna ci
g�w jest wielko[
ci
definiowan
jako: Rxy(m) = E{x(n) �" y(n + m)} ( 0 ) gdzie x(n) i y(n) s
stacjonarnymi procesami losowymi, - " <� n <� " oraz E jest operatorem warto[
ci oczekiwanej. Kowariancja wzajemna jest r�wna korelacji z usuni
t
[
redni
. Cxy(m) = E�� (x(n) - mx )(y(n + m) - m )�� �� �� ( 1 ) y �� �� gdzie m i m s
warto[
ciami [
rednimi korelowanych ci
g�w. x y ZwykB
a estymata bazuj
ca na M pr�bkach x(n) i y(n) jest deterministyczn
korelacj
wzajemn
ci
g�w: M- m - 1 �� x(n)y(n + m) dla m e" 0 �� " � Rxy(m) = �� n= 0 ( 2 ) �� � Rxy(- m) dla m <� 0 �� ZakB
ada si
, |
e x(n) i y(n) s
indeksowane od 0 do M - 1, a R (m), od -(M - 1) do M - 1. xy Funkcja xcorr realizuje powy|
szy wz�r, wykorzystuj
c algorytmy bazuj
ce na FFT dla danych ci
g�w wej[
ciowych x(n) i y(n), zawartych w wektorach x i y o dB
ugo[
ciach M. PrzykB
adowo: x=[1 1 1 1 1]'; y=x; xyc=xcorr(x,y) " xyc = " 1.0000 " 2.0000 " 3.0000 " 4.0000 " 5.0000 " 4.0000 " 3.0000 " 2.0000 " 1.0000 Sekwencja wynikowa ma dB
ugo[

dwa razy wi
ksz
od dB
ugo[
ci ci
gu wej[
ciowego minus jeden. Tak wi
c M-ty element jest korelowany przy zerowej zwB
oce. R�wnie|
zauwa|
amy znane impulsy tr�jk
tne na wyj[
ciu, b
d
ce wynikiem splotu dw�ch przebieg�w kwadratowych. Funkcja xcov estymuje kowariancj
wB
asn
i wzajemn
danych wej[
ciowych. Posiada ona te same opcje i przelicza te same r�wnania, co funkcja korelacji, ale najpierw usuwa warto[
ci [
rednie z przebieg�w x i y. Obci
|
enie i normalizacja Estymator jest obci
|
ony, je[
li jego warto[

oczekiwana nie jest r�wna warto[
ci wielko[
ci, kt�r
estymuje. Warto[

oczekiwana dla przebiegu funkcji xcorr wynosi: M- m - 1 )� E(R (m)) = E{x(n)y(n + m)} = (M - m )R (m) ( 3 ) xy " xy n= 0 Funkcja xcorr wyznacza estymator nieobci
|
ony, dziel
c sum
(3) przez M-��m��, je[
li zadamy parametr  unbiased po sekwencji element�w wej[
ciowych: xcorr(x,y, unbiased ); Mimo |
e estymator jest nieobci
|
ony, punkty koD
cowe (bliskie -(M-1) i M-1 maj
stosunkowo du|

wariancj
, poniewa|
funkcja xcorr wylicza je na podstawie niewielu punkt�w danych. Mo|
liwym rozwi
zaniem jest dzielenie estymaty przez M, u|
ywaj
c w tym celu parametru  biased : xcorr(x,y, biased ); W tym schemacie jedynie pr�bka korelacji przy zerowym op�z
nieniu (M-ty element wyj[
cia) jest nieobci
|
ona. Ta estymata jest cz
sto bardziej po|

dana ni|
nieobci
|
ona, poniewa|
unika ona losowej, du|
ej wariancji w koD
cowych punktach korelowanej sekwencji. Funkcja xcorr umo|
liwia tak|
e zastosowanie innego schematu normalizacji: xcorr(x,y,'coeff'); Polega on na podzieleniu estymaty wyj[
ciowej przeznorm(x)norm(y) w taki spos�b, |
e wsp�B
czynnik korelacji przy zerowym op�z
nieniu ma warto[

jeden. SygnaB
y wielokanaB
owe Dla sygnaB
u wielokanaB
owego funkcje xcorr i xcov estymuj
wB
asn
i wzajemn
korelacj
i kowariancj
ci
gu danych dla wszystkich kanaB
�w jednocze[
nie. Maj
c dan
macierz S o wymiarach M na N, reprezentuj
c
sygnaB
N kanaB
owy (jedna kolumna dla pojedynczego kanaB
u), funkcja xcorr(S) oblicza macierz korelacji wB
asnej i wzajemnej dla sygnaB
�w zawartych w S. Otrzymana macierz wynikowa posiada wymiary (2M-1) na N2. PrzykB
adowo, je[
li S jest sygnaB
em trzykanaB
owym S=[s1 s2 s3]; to rezultaty dziaB
ania funkcji R=xcorr(S) s
podawane nast
puj
co: R=[Rs1s1 Rs1s2 Rs1s3 Rs2s1 Rs2s2 Rs2s3 Rs3s1 Rs3s2 Rs3s3]. W Matlabie dost
pne s
ponadto dwie funkcje cov i corrcoef. Stosuj
c te funkcje, mo|
na wyznaczy
estymaty kowariancji oraz znormalizowanej kowariancji pomi
dzy r�|
nymi kanaB
ami przy zerowym przesuni
ciu. Wy|
ej wymienione funkcje zestawiaj
wyniki swego dziaB
ania w macierzy kwadratowej. G
sto[

widmowa mocy Celem analizy widmowej jest wyznaczenie cz
stotliwo[
ciowej zawarto[
ci sygnaB
u, procesu losowego lub ukB
adu, bazuj
c na skoD
czonych zbiorach danych. Estymacja widmowej g
sto[
ci mocy jest u|
yteczna w r�|
nych zastosowaniach, wB

czaj
c detekcj
sygnaB
�w zakB
�conych szumem o szerokim pa[
mie. Widmowa g
sto[

mocy (PSD) stacjonarnego procesu losowego x(n) jest zwi
zana z funkcj
korelacji sygnaB
u dyskretnego poprzez dyskretn
transformacj
Fouriera. " Pxx (� ) = R (m)e- j� m ( 4 ) " xx m= - " Funkcja g
sto[
ci widmowej mocy posiada tak
wB
asno[

, |
e jej caB
ka w danym pa[
mie cz
stotliwo[
ci jest r�wna mocy sygnaB
u x(n) w tym pa[
mie. PSD jest specjalnym przypadkiem funkcji wzajemnej g
sto[
ci widmowej (CSD), zdefiniowanej pomi
dzy dwoma sygnaB
ami x(n) i y(n) w postaci wzoru: " Pxy (� ) = R (m)e- j� m ( 5 ) " xy m= - " Tak jak w wypadku funkcji korelacji i kowariancji, polecenia Matlaba estymuj
PSD i CSD dla sygnaB
�w o skoD
czonej dB
ugo[
ci. Istnieje wiele r�|
nych metod estymacji PSD. Mo|
na jednak zgrupowa
je w dw�ch kategoriach: estymacji parametrycznej i nieparametrycznej. GB
�wn
technik
stosowan
w SPT jest nieparametryczny algorytm, opracowany przez Welcha. W przyborniku wykorzystano t
metod
w funkcjach psd, csd, tfe, oraz cohere. W[
r�d metod parametrycznych stosowana jest funkcja pmem, dokonuj
ca estymacji widma metod
 maksymalnej entropii . Bardziej szczeg�B
owe om�wienie funkcji do estymacji PSD, wykorzystuj
cych algorytmy parametryczne, znajduje si
plikach pomocy Matlaba. Metoda Welcha Jedn
z metod estymacji g
sto[
ci widmowej mocy jest wyznaczenie dyskretnej transformaty Fouriera pr�bek konkretnej realizacji danego procesu, a nast
pnie obliczenie kwadratu z warto[
ci bezwzgl
dnej wyniku. PrzykB
adowo, wyznaczono PSD dla 1001-elementowego sygnaB
u xn, skB
adaj
cego si
z dw�ch sinusoid i szumu. Fs=1000; t=0:1/Fs:1; xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+... randn(size(t)); Zgrubny estymator PSD wyznacza si
z xn stosuj
c instrukcje: Pxx=abs(fft(xn, 1024)).^2/1001; Powy|
szy estymator jest nazywany periodogramem. Skalowanie kwadratu warto[
ci bezwzgl
dnej z FFT przez kwadrat normy okna danych zastosowanego do sygnaB
u (w tym wypadku okna prostok
tnego o dB
ugo[
ci 1001), gwarantuje to, |
e estymator ten jest asymptotycznie nieobci
|
ony. To znaczy, |
e je[
li liczba pr�bek wzrasta, warto[

oczekiwana periodogramu przybli|
a si
do prawdziwej warto[
ci PSD. Problemem przy estymacji PSD metod
periodogramu jest fakt, |
e jego wariancja jest du|
a i nie maleje przy wzrastaj
cej liczbie pr�bek. Nast
pne dwa przykB
ady pokazuj
, |
e pomimo wzrostu dB
ugo[
ci analizowanego wektora, periodogram nie wygB
adza si
. Pxx_short=abs(fft(xn,256)).^2/256; subplot(2,1,1) plot((0:255)/256*Fs,10*log10(Pxx_short),'black') grid, title('Periodogram kr�tkiego sygnaB
u') xlabel('Cz
stotliwo[

w Hz'), ylabel('Widmo mocy (dB)'); subplot(2,1,2), plot((0:1023)/1024*Fs,10*log10(Pxx),'black') grid, text(250,45,'Periodogram dB
u|
szego sygnaB
u') xlabel('Cz
stotliwo[

w Hz'), ylabel('Widmo mocy (dB)'); Periodogram kr�tkiego sygnaB
u 40 20 0 -20 -40 0 200 400 600 800 1000 Cz
stotliwo[

w Hz Periodogram dB
u|
szego sygnaB
u 40 20 0 -20 -40 0 200 400 600 800 1000 Cz
stotliwo[

w Hz Rys. 1 Periodogramy dla sekwencji o r�|
nych dB
ugo[
ciach Redukcja wariancji estymatora PSD nast
puje przez podzielenie sygnaB
u na nie nakB
adaj
ce si
fragmenty oraz u[
rednianie ich widm w dziedzinie cz
stotliwo[
ci. subplot(1,1,1), Pxx=(abs(fft(xn(1:256))).^2+... abs(fft(xn(257:512))).^2+... abs(fft(xn(513:768))).^2)/(256*3); plot((0:255)/256*Fs,10*log10(Pxx),'black') grid, xlabel('Cz
stotliwo[

w Hz'), ylabel('Widmo mocy (dB)') title('U[
redniony periodogram (fragmenty sygnaB
u nie nakB
adaj
si
)') U[
redniony periodogram (fragmenty sygnaB
u nie nakB
adaj
si
) 25 20 15 10 5 0 -5 -10 0 200 400 600 800 1000 Cz
stotliwo[

w Hz Rys. 2 Periodogram u[
redniony w nie nakB
adaj
cych si
przedziaB
ach Taki u[
redniony estymator posiada jedn
trzeci
wariancji periodogramu o dB
ugo[
ci 256, otrzymanego wcze[
niej. U[
redniaj
c wi
cej pr�bek, mo|
na otrzyma
w rezultacie ni|
sz
wariancj
. DB
ugo[

sygnaB
u ogranicza jednak|
e liczb
mo|
liwych do utworzenia odcink�w Widmo mocy (dB) Widmo mocy (dB) Widmo mocy (dB) (do 3 sekcji o dB
ugo[
ci 256 w poprzednim przykB
adzie). By uzyska
wi
ksz
ilo[

sekcji, dzielimy sygnaB
na fragmenty nakB
adaj
ce si
. Pxx=(abs(fft(xn(1:256))).^2+... abs(fft(xn(129:384))).^2+... abs(fft(xn(257:512))).^2+... abs(fft(xn(385:640))).^2+... abs(fft(xn(513:768))).^2+... abs(fft(xn(641:896))).^2)/(256*6); plot((0:255)/256*Fs,10*log10(Pxx),'black'),grid, xlabel('Cz
stotliwo[

w Hz'), ylabel('Widmo mocy (dB)') title('U[
redniony periodogram (pr�bki nakB
adaj
si
)') W tym wypadku odcinki s
statystycznie zale|
ne, co z kolei powoduje wzrost wariancji. Nale|
y wi
c znalez

kompromis pomi
dzy liczb
fragment�w sygnaB
u i stopniem ich nakB
adania. U[
redniony periodogram (pr�bki nakB
adaj
25 si
) 20 15 10 5 0 -5 -10 0 200 400 600 800 1000 Cz
stotliwo[

w Hz Rys. 3 Periodogram u[
redniony w nakB
adaj
cych si
przedziaB
ach Inn
drog
do poprawienia wB
asno[
ci statystycznych estymatora periodogramu jest zastosowanie nieprostok
tnego okna czasowego do wygB
adzenia koD
c�w analizowanych danych. W rezultacie otrzymujemy zmodyfikowany periodogram. Operacja ta redukuje efekt zale|
no[
ci nakB
adaj
cych si
sekcji, poniewa|
okno zmniejsza amplitud
do zera przy ich brzegach. Nieprostok
tne okno czasowe zmniejsza poziom bocznych interferencji (lub  przecieku widma ), r�wnocze[
nie jednak powoduj
c wzrost szeroko[
ci pr
|
k�w widmowych. Z odpowiednim oknem czasowym (Hamminga, Hanninga lub Kaisera) oraz 50 procentowym nakB
adaniem si
analizowanych odcink�w sygnaB
u uzyskuje si
estymator o znacznie zmniejszonej wariancji. Rezultat zastosowania okna Hanninga do analizowanego sygnaB
u przedstawiono na rys. 4. W tym celu zrealizowano nast
puj
ce funkcje: w=hanning(256)'; Pxx=(abs(fft(w.*xn(1:256))).^2+... abs(fft(w.*xn(129:384))).^2+... abs(fft(w.*xn(257:512))).^2+... abs(fft(w.*xn(385:640))).^2+... abs(fft(w.*xn(513:768))).^2+... abs(fft(w.*xn(641:896))).^2)/(norm(w)^2*6); plot((0:255)/256*Fs,10*log10(Pxx),'black') Widmo mocy (dB) grid, xlabel('Cz
stotliwo[

w Hz'), ylabel('Widmo mocy (dB)') title( Periodogram zmodyfikowany oknem Hanninga') Periodogram zmodyfikowany oknem Hanninga 25 20 15 10 5 0 -5 -10 0 200 400 600 800 1000 Cz
stotliwo[

w Hz Rys. 4 Periodogram zmodyfikowany oknem Hanninga Na rysunku 4 mo|
na zauwa|
y
, |
e pr
|
ki widma rozszerzyB
y si
, a poziom szum�w wydaje si
by
mniejszy w por�wnaniu z poprzednimi wykresami. Przedstawione tutaj metody: u[
redniania oraz zmodyfikowanego periodogramu nale|

do metod opartych na algorytmie Welcha. Zawarte w przyborniku SPT funkcje psd i csd umo|
liwiaj
dob�r odpowiednich parametr�w analizy (dB
ugo[

FFT, rodzaj okna, szeroko[

nakB
adania si
pr�bek) w trakcie obliczeD
PSD i CSD dla jednego lub dw�ch sygnaB
�w przy pomocy metody Welcha. Obliczanie funkcji widmowej g
sto[
ci mocy za pomoc
procedur SPT Funkcja psd u[
rednia i skaluje zmodyfikowany periodogram z realizacji sygnaB
u pozbawionych trendu. Jako argumenty funkcji nale|
y wyszczeg�lni
parametry, kt�re steruj
wykonaniem algorytmu. Funkcja psd, obliczaj
c estymator PSD dla sekwencji xn u|
ywa domy[
lnie okna czasowego Hanninga oraz FFT o dB
ugo[

256 pr�bek. Zadany ci
g jest dzielony na nienakB
adaj
ce si
sekcje oraz jest usuwany trend. Pxx=psd(xn); Je[
li oryginalna sekwencja wyra|
ona jest w woltach V, Pxx posiada jednostk
V2. Poni|
ej zostaB
powt�rzony ostatni z podanych przykB
ad�w. Do jego realizacji przyj
to 128 jako liczb
nakB
adaj
cych si
pr�bek oraz nie usuni
to trendu z sygnaB
u. nfft=256; Fs=1000; window=hanning(256); noverlap=128; dflag='none'; Pxx=psd(xn,nfft,Fs,window,noverlap,dflag); Cz
stotliwo[

pr�bkowania zadawana jako jeden z parametr�w, nie odgrywa wi
kszej roli w obliczeniach, pomaga jednak w skalowaniu osi cz
stotliwo[
ci wykresu. Funkcja psd bez |
adnych parametr�w wyj[
ciowych generuje wykres przedstawiony na rys 5. psd(xn,nfft,Fs,window,noverlap,dflag); xlabel('Cz
stotliwo[

Hz'), ylabel('Amplituda widmowa mocy (dB)') Widmo mocy (dB) 25 20 15 10 5 0 -5 -10 0 100 200 300 400 500 Cz
stotliwo[

Hz Rys. 5 Periodogram uzyskany w wyniku dziaB
ania funkcji psd Dost
p do wektora, dla kt�rego obliczana jest PSD, uzyskujemy dodaj
c drugi parametr wyj[
ciowy. Wektor ten jest u|
yteczny podczas nieautomatycznego wykonywania wykresu PSD. [Pxx,f]=psd(xn, Fs, 256, noverlap, dflag); plot(f,10*log10(Pxx)) Dla rzeczywistego sygnaB
u xn, psd oblicza widmowe g
sto[
ci mocy jedynie dla cz
stotliwo[
ci od 0 do cz
stotliwo[
ci Nyquista. We wcze[
niejszych przykB
adach, w kt�rych estymator PSD byB
generowany przy pomocy FFT, zakres cz
stotliwo[
ci zawieraB
si
w granicach od 0 do cz
stotliwo[
ci pr�bkowania. Obci
|
enie i normalizacja Analizuj
c przedstawione powy|
ej wyniki obliczeD
funkcji psd, mo|
na zauwa|
y
kilka charakterystycznych cech analizowanego sygnaB
u xn. Poziom szum�w jest maB
y, oscyluje wok�B
zera decybeli (dB) i z zaB
o|
enia jest to biaB
y szum o wariancji 1. Ponadto moc sygnaB
u xn jest skoncentrowana w dw�ch pr
|
kach o cz
stotliwo[
ciach 50 i 120 Hz. Pr
|
ek o cz
stotliwo[
ci 50 Hz jest 6 dB poni|
ej pr
|
ka o cz
stotliwo[
ci 120 Hz. Na tej podstawie mo|
na stwierdzi
, |
e sinusoida o wy|
szej cz
stotliwo[
ci posiada dwukrotnie wi
ksz
amplitud
(106/20=2.0). Niestety na podstawie rzeczywistej wysoko[
ci pr
|
k�w nie mo|
na wiele powiedzie
o oryginalnej amplitudzie sinusoid bez dokonania dodatkowych analiz. W celu uzyskania u|
ytecznych informacji o amplitudzie pr
|
k�w badanych sinusoid, trzeba zauwa|
y
, |
e warto[

oczekiwana estymatora PSD jest okre[
lana wzorem: � )� 1 2 E{Pxx (� )} = ( 6 ) xx 2 +"P (� ) W(� - � ) d� 2� w - � Wielko[

ta, jest splotem warto[
ci PSD z podniesion
do kwadratu warto[
ci
bezwzgl
dn
dyskretnej transformaty Fouriera zastosowanego okna W(�), podzielony przez kwadrat jego normy. Wielko[

skaluj
ca (norma okna) jest sum
kwadrat�w wsp�B
czynnik�w funkcji okna: 2 w = w(n)2 ( 7 ) " R�wnanie (6) wskazuje na fakt, |
e je[
li P (�) posiada pr
|
ek o wysoko[
ci 1 przy okre[
lonej xx 2 2 W(0) / w cz
stotliwo[
ci � , to estymator b
dzie w przybli|
eniu przyjmowaB
warto[

dla 0 tej cz
stotliwo[
ci. Tak wi
c, by uzyska
estymator amplitudy oryginalnych pr
|
k�w, nale|
y pomno|
y
wynik dziaB
ania funkcji psd przez norm2(w)/sum2(w), gdzie w jest wektorem okna. Amplituda widmowa mocy (dB) Wynika to z twierdzenia Parsevala i|
sum2(w)=sum2(W). Skalowanie jest niezale|
ne od dB
ugo[
ci i ksztaB
tu okna. PrzykB
adowo: w1=hanning(256); w2=hanning(500); [Pxx1,f1]=psd(xn,256,Fs,w1,128,'none'); [Pxx2,f2]=psd(xn,1024,Fs,w2,250,'none'); subplot(2,1,1), plot(f1,10*log10(Pxx1*norm(w1)^2/sum(w1)^2),... 'black'),grid, subplot(2,1,2), plot(f2,10*log10(Pxx2*norm(w2)^2/sum(w2)^2),... 'black'),grid Na obydwu wykresach (rys. 6), kt�re pokazuj
widma jedynie dla dodatnich cz
stotliwo[
ci (cz
stotliwo[
ci ujemne s
takie same), pr
|
ek dla wy|
szej cz
stotliwo[
ci posiada warto[

0 dB, a dla ni|
szej -6 dB. Amplituda sinusoidy o cz
stotliwo[
ci 120 Hz, 0 dB, odpowiada kwadratowi amplitudy r�wnemu 1. Wynika to z faktu, i|
sinusoida o amplitudzie r�wnej 2 posiada widmo z dwoma eksponencjalnymi skB
adnikami o amplitudzie 1, dla dodatnich i ujemnych cz
stotliwo[
ci. Podobnie, sinusoida o cz
stotliwo[
ci 50 Hz posiada dwa skB
adniki cz
stotliwo[
ci, dodatni i ujemny, r�wne kwadratom amplitudy (1/2)2 lub 10log10(0.25)= -6 dB. Nale|
y r�wnie|
zauwa|
y
, |
e drugi wykres posiada nieznacznie ni|
szy poziom szumu, b
d
cy wynikiem wi
kszej dB
ugo[
ci okna. 0 - 1 0 - 2 0 - 3 0 - 4 0 0 1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 0 - 1 0 - 2 0 - 3 0 - 4 0 0 1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 Rys. 6 Przeskalowane periodogramy Zwi
zek Parsevala CaB
ka z PSD, obliczona dla caB
ego pasma cz
stotliwo[
ci, jest miar
caB
kowitej energii sygnaB
u. Za pomoc
obliczenia funkcji psd mo|
na w przybli|
eniu sprawdzi
t
zale|
no[

(zgodno[

obliczeD
energii w dziedzinie czasu i cz
stotliwo[
ci) dla danego sygnaB
u: [Pxx,f]=psd(xn,256,Fs,256,128,'none'); format long sum(xn.^2)/length(xn) " ans = " 3.60376766505040 sum(Pxx)/length(Pxx) " ans = " 3.68430077838700 Obliczenie procentowej zawarto[
ci energii w sygnale w danym pa[
mie cz
stotliwo[
ci, wymaga zsumowania estymatora PSD jedynie dla wybranych pr
|
k�w widma. Tak wi
c procentow
energi
sygnaB
u xn w pa[
mie od 40 do 60 Hz wyznacza si
stosuj
c sekwencje instrukcji Matlaba: ind=find(f>40&f<�60) " ind = " 12 " 13 " 14 " 15 " 16 sum(Pxx(ind))/sum(Pxx) " ans = " 0.13369450463404 Obliczenie funkcji wzajemnej g
sto[
ci widmowej Estymuj
c, przy pomocy metody Welcha, wzajemn
g
sto[

widmow
dw�ch sygnaB
�w o r�wnej dB
ugo[
ci x i y, funkcja csd oblicza periodogram jako iloczyn FFT z x oraz warto[
ci sprz
|
onej do wyniku FFT z y. Przetwarzanie sygnaB
�w przy pomocy funkcji csd (a wi
c: dzielenie na pr�bki, usuwanie trendu, stosowanie okien) odbywa si
w ten sam spos�b, jak dla funkcji psd. Pxy=csd(x,y,nfft,window,nowerlap,dflag); PrzedziaB
y ufno[
ci Za pomoc
funkcji psd i csd mo|
na oblicza
przedziaB
y ufno[
ci dla estymator�w g
sto[
ci widmowych. Mo|
na wprowadzi
jako argument wej[
ciowy parametr p, kt�ry podaje przedziaB
ufno[
ci (w procentach). [Pxx,Pxxc,f]=psd(x,nfft,Fs,window,noverlap,... p,dflag); [Pxy,Pxyc,f]=csd(x,y,nfft,Fs,window,noverlap,... p,dflag); Parametr p musi by
wielko[
ci
skalarn
zawart
w granicach od 0 do 1. Funkcje zakB
adaj
normalny rozkB
ad periodogramu z nienakB
adaj
cymi si
pr�bkami. Je[
li te zaB
o|
enia s
speB
nione, wtedy p�"100% jest prawdopodobieD
stwem, |
e w przedziale[Pxx- Pxxc(:,1), Pxx+Pxxc(:,2)] zawiera si
rzeczywista warto[

PSD. Je[
li pr�bki nakB
adaj
si
, przedziaB
ufno[
ci nie jest wyznaczony poprawnie, a funkcja wy[
wietla ostrze|
enie. Innym rozwi
zaniem jest przyj
cie zaB
o|
enia, |
e periodogramy posiadaj
rozkB
ad chi kwadrat. Je[
li dost
pny jest przybornik Statistics Toolbox, mo|
na obliczy
poziom ufno[
ci p�"100% dla tak przyj
tego rozkB
adu, stosuj
c nast
puj
ce funkcje: k=fix((length(x)-noverlap)/(length(window)-noverlap)); alfa=1-p; confid=2*k*Pxx*(1./chi2inv([1-alfa/2 alfa/2],2*k)); gdzie k jest liczb
sekcji u|
ytych do obliczenia estymatora PSD, za[
confid jest dwukolumnow
macierz
, skB
adaj
c
si
z g�rnych i dolnych zakres�w ufno[
ci. W�wczas mamy p�"100%-owe prawdopodobieD
stwo, |
e poziom ufno[
ci [Pxx - confid(:,1) Pxx+confid(:,2)] zawiera rzeczywist
PSD. Ten poziom ufno[
ci jest wiarygodny jedynie, gdy sekcje nie pokrywaj
si
. Estymator funkcji przej[
cia Jednym z zastosowaD
metody Welcha jest identyfikacja modeli funkcjonalnych. Jednym z tego rodzaju modeli jest funkcja przej[
cia ukB
adu. ZakB
ada si
, |
e badany obiekt H jest ukB
adem liniowym o staB
ych wsp�B
czynnikach, oraz x(n) i y(n) s
odpowiednio jego wej[
ciem i wyj[
ciem. Wtedy prawdziwa jest zale|
no[

: Pxy (� ) = H(� )Pxx (� ) ( 8 ) gdzie H jest widmow
funkcj
przej[
cia. Znaj
c sygnaB
y x(n) i y(n) mo|
na zdefiniowa
estymator funkcji przej[
cia ukB
adu jako: � Pxy (� ) $
(� ) = ( 9 ) � Pxx (� ) Metoda ta pozwala na estymacj
zar�wno amplitudy jak i fazy sygnaB
u. Funkcja fte wykorzystuje metod
Welcha do obliczenia CSD i PSD, a nast
pnie formuje ich iloraz w celu estymowania funkcji przej[
cia. Funkcja tfe posiada tak
sam
skB
adni
jak funkcja csd. W poni|
szym programie przefiltrowano sygnaB
xn przy pomocy filtru FIR, a nast
pnie wyznaczono charakterystyk
amplitudowo-cz
stotliwo[
ciow
oraz estymowan
funkcj
przej[
cia: h=ones(1,10)/10; yn=filter(h,1,xn); [Hest,f]=tfe(xn,yn,256,Fs,256,128,'none'); H=freqz(h,1,f,Fs); subplot(2,1,1), plot(f,abs(H),'black') title('Rzeczywista amplituda funkcji przej[
cia') xlabel('Cz
stotliwo[

w Hz') subplot(2,1,2), plot(f,abs(Hest),'black') text(110,1.08,'Estymata amplitudy funkcji przej[
cia') xlabel('Cz
stotliwo[

w Hz') Funkcja koherencji Funkcja koherencji okre[
la zwi
zek pomi
dzy dwoma sygnaB
ami w dziedzinie cz
stotliwo[
ci. Kwadrat moduB
u koherencji pomi
dzy sygnaB
ami x(n) i y(n) wyznacza si
ze wzoru: 2 Pxy (� ) ( 10 ) Cxy (� ) = Pxx (� )Pyy (� ) Iloraz ten jest liczb
rzeczywist
z przedziaB
y [0,1] i jest miar
korelacji pomi
dzy sygnaB
ami x(n) i y(n) dla danej cz
stotliwo[
ci �. Funkcja cohere dla zadanych ci
g�w pr�bek x i y oblicza PSD i CSD, a nast
pnie iloraz kwadratu warto[
ci bezwzgl
dnej CSD przez iloczyn PSD dla x i y. Jej parametry i spos�b dziaB
ania jest podobny do funkcji csd i tfe. Funkcj
koherencji pomi
dzy ci
giem wej[
ciowym xn i wyj[
ciowym yn filtru rozwa|
anego w poprzednim rozdziale wyznaczono za pomoc
sekwencji instrukcji: Rzeczywista amplituda funkcji przej[
cia 1 0.5 0 0 100 200 300 400 500 Cz
stotliwo[

w Hz Estymata amplitudy funkcji przej[
cia 1 0.5 0 0 100 200 300 400 500 Cz
stotliwo[

w Hz Rys. 7 Charakterystyki amplitudowo-cz
stotliwo[
ciowe zadanej i estymowanej funkcji przej[
cia filtru clf, cohere(xn,yn,256,Fs,256,128,'none') title('Funkcja koherencji') Funkcja koherencji 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 100 200 300 400 500 Cz
stotliwo[

Rys. 8 Przebieg funkcji koherencji Dla sygnaB
�w z poprzedniego przykB
adu koherencja obliczona za pomoc
funkcji cohere wynosi 1. Uwaga W obecnych wersjach toolboxu Signal Processing przepisano funkcje obliczaj
cych estymatory dziaB
aj
ce w dziedzinie cz
stotliwo[
ci. W chwili obecnej zmiennie mog
by
stosowane funkcje: Estymator funkcji koherencji  do obliczania g
sto[
ci widmowej mocy:  [Pxx,f] = psd(xn,nfft,Fs,window,noverlap,dflag);  [Pxx,f] = pwelch(xn,window,noverlap,nfft,Fs);  do obliczania wzajemnej g
sto[
ci widmowej  [Pxx,f] = csd(xn,yn,nfft,Fs,window,noverlap,dflag)  [Pxx,f] = cpsd(xn,yn,window,noverlap,nfft,Fs);  do obliczania widmowej funkcji przej[
cia:  [Pxx,f] = tfe(xn,yn,nfft,Fs,window,noverlap,dflag)  [Pxx,f] = tfestimate(xn,yn,window,noverlap,nfft,Fs);  do obliczania funkcji koherencji:  [Pxx,f] = cohere(xn,yn,nfft,Fs,window,noverlap,dflag)  [Pxx,f] = mscohere(xn,yn,window,noverlap,nfft,Fs); Zadania do wykonania 1. Wykona
symulacj
ukB
adu testowanego o dw�ch stopniach swobody dla parametr�w ustalonych w trakcie zaj

z modelowania. Zarejestruj sygnaB
y testowe w postaci przebieg�w wymuszaj
cych oraz sygnaB
�w odpowiedzi dla poszczeg�lnych stopni swobody. SygnaB
y odpowiedzi powinny by
w postaci wektor�w przemieszczeD
, pr
dko[
ci i przyspieszeD
o odpowiedniej dB
ugo[
ci. Dane: M = 10, M = 7, C = 0.7, C = 0.9, C = 0.7, K = 120, K = 200, K = 170 1 2 1 2 3 1 2 3 2. Dla zapisanych danych dokonaj wyznaczenia estymator�w g
sto[
ci widmowych mocy, widmowych funkcji przej[
cia oraz funkcji koherencji. Aby estymatory miaB
y znaczenie statystyczne konieczna jest praca z przebiegami o odpowiedniej dB
ugo[
ci. Wykonaj niezale|
nego wyznaczenia stosownych estymator�w bazuj
c na sygnale wymuszaj
cym oraz: " sygnaB
ach przyspieszeD
" sygnaB
ach pr
dko[
ci " sygnaB
ach przemieszczeD
. 3. Por�wnaj wykonuj
c stosowne wykresy wyniki estymacji uzyskane w punkcie 2. Skomentuj wyniki eksperymentu. 4. Powt�rz zadania 2 i 3 dodaj
c do symulacyjnych danych wyj[
ciowych szum o amplitudach na poziomie 10% i 20% sygnaB
u.