Podobieństwo sygnałów – korelacja
Iloczyn skalarny wektorów/sygnałów
W przestrzeni L
2
ℝ
〈
x , y 〉=
∫
x t y t dt
W przestrzeni ℝ
N
〈
x , y 〉=
1
N
∑
n=0
N −1
x n y n
Jeżeli 〈 x , y 〉=0 to x ⊥ y
a co jeżeli
〈
x , y 〉≠0
lub inaczej
∣〈
x , y 〉∣0
?
Korelacja
Korelacja w przestrzeni
L
2
ℝ
w przypadku stacjonarnym
R=
∫
x t x t− dt
po dyskretyzacji w przestrzeni ℝ
N
Rl =lim
n ∞
1
N
∑
n=0
N −1
x n x n−l
Korelacja w przestrzeni
ℝ
N
Z teorii procesów stochastycznych
R
x
k , m=
1
2
E[ x k − x m−]
gdzie E [.] oznacza operator wartości oczekiwanej (w dużym uproszczeniu jest to wartość
średnia)
wartość średnia procesu losowego
wariancja procesu losowego
Konsekwentnie
R
xy
k , m=
1
x
y
E [ x k −
x
y m−
y
]
Zwykle zakładamy że:
●
=
0
●
=
1
●
proces losowy (sygnał) jest stacjonarny wtedy
x n , x n−l
lub
x n , y n−l
- 1 -
Użyteczne definicje
R
x
l= E [x n x n−l] - autokorelacja
R
xy
l =E [ x n y n−l]
- korelacja wzajemna (kroskorelacja)
W praktyce można różnie liczyć estymator wartości oczekiwanej
R
x
l=
1
N
∑
n
x n x n−l - estymator obciążony
R
x
l=
1
N −l
∑
n=0
N −l −1
x n x n−l - estymator nieobciążony
Wyjaśnić pojęcia:
●
współczynnik korelacji
●
unormowany współczynnik korelacji ( 1/
x
2
)
●
funkcja korelacji
●
unormowana funkcja korelacji
●
miara podobieństwa sygnałów (dla l=0 otrzymujemy iloczyn skalarny !!!)
Przykład:
N=1000;n=(0:N-1);x=sin(2*pi*5/N*n+.3*pi)+randn(1,N);plot(n,x);
s1=sin(2*pi*5/N*n);s2=sin(2*pi*4/N*n);s3=sin(2*pi*13/N*n);s4=sin(2*pi*10/N*n)
;s5=sin(2*pi*4.8/N*n);
max(abs(xcorr(x,s1)))
max(abs(xcorr(x,s2)))
max(abs(xcorr(x,s3)))
max(abs(xcorr(x,s4)))
max(abs(xcorr(x,s5)))
Własności funkcji autokorelacji
1.
R
x
l= R
x
−
l
funkcja parzysta
2.
R
x
0≥R
x
l wartość maksymalna dla zerowego przesunięcia
- 2 -
Transformacje czasowo-częstotliwościowe
Krótkoterminowa transformata Fouriera (ang. STFT)
STFT {x n}≡ X k , l =
∑
n =0
N −1
x n w n−l e
−
j 2
N
k n
gdzie l – dyskretny czas, k – dyskretna częstotliwość
spektrogram to:
S k , l=∣X k , l∣
2
Przykład:
N=64;n=(0:N-1);
x=1.2*sin(2*pi*.13*n);y=2*sin(2*pi*.07*n);z=.8*sin(2*pi*.27*n);s=[x,y,z];
s=s+.1*randn(size(s));plot(s);
S=fft(s); M=size(S,2); f=(0:M-1)./M; plot(f,abs(S));
w=gausswin(N,3)';plot(w);
okno1=[w,zeros(1,N),zeros(1,N)];
okno2=[zeros(1,N),w,zeros(1,N)];
okno3=[zeros(1,N),zeros(1,N),w];
plot(okno1);hold on;plot(okno2);plot(okno3);hold off;
plot(s.*okno1);hold on;plot(s.*okno2);plot(s.*okno3);hold off;
S1=fft(s.*okno1); S2=fft(s.*okno2); S3=fft(s.*okno3);
M=size(S1,2); f=(0:M-1)./M;
plot(f,abs(S1),f,abs(S2),f,abs(S3));
N=256;n=(0:N-1);
x=sin(2*pi*.13*n);y=sin(2*pi*.07*n);z=sin(2*pi*.27*n);s=[x,y,z];
s=s+.7*randn(size(s));plot(s);
M=128; w=gausswin(M,3)'; plot(w);
tmp=[zeros(1,M/2),s,zeros(1,M/2)];plot(tmp);
L=3*N;S=zeros(L,M);
for l=(1:L), v=tmp(l:l+M-1).*w;V = fft(v);S(l,:)=abs(V).^2;end;
l=(0:L-1);f=(0:M-1)./M;imagesc(l,f,S');
Transformata Wignera-Villa
X
WV
k , l=
∑
n=0
N −1
x l x l−ne
−
j 2
N
k n
!!! Problemy z częstotliwością Nquista !!!
N=100;n=(0:N-1);x=sin(2*pi*.13*n);y=zeros(1,N);z=sin(2*pi*.28*n);s=[x,y,z];
plot(s);
L= 3*N; l=(0:L-1); M=256; [TFR,nt,nf]=tfrwv(s',n,M);imagesc(nt,nf,TFR);
N=500;n=(0:N-1);f=linspace(.1,.3,N);x=sin(2*pi*f.*n);
M=128;[TFR,nt,nf]=tfrwv(x',n,M);imagesc(nt,nf,TFR);
- 3 -
M=128;[TFR,nt,nf]=tfrwv(hilbert(x'),n,M);imagesc(nt,nf,TFR);
Dyskretna transformata falkowa (ang. wavelet)
Wavelet/Falka – (mała fala) sygnał okresowy szybko zanikający do zera
Falka ciągła:
a , b
t=
1
a
t−b
a
gdzie
t ∈L
2
ℝ
,
a , b∈ℝ ,
a0
t prototyp falki, falka matka
–
b przesunięcie w czasie
–
a skalowanie w częstotliwości
Stąd nazwa transformacja „częstotliwość-skala” (skalogram)
Transformata falkowa to iloczyn skalarny badanego sygnału
z prototypami falek
Reprezentacja w przestrzeni
L
2
ℝ
W a ,b =
〈
x t ,
a , b
t
〉
=
∫
x t
a ,b
t dt=
1
a
∫
ℝ
x t
t−b
a
dt
przesunięcie w czasie:
y t= x t−u ,
W
y
a ,b=W
x
a ,b−u
przesunięcie w częstotliwości
y t= x st ,
W
y
a , b=
1
s
W
x
sa , sb
Transformata odwrotna:
x t =
∫
ℝ
∫
ℝ
W a , b
a ,b
t da db
Mamy bazę ortonormalną – rodzina falek, mamy współczynniki reprezentacji, mamy też
transformatę odwrotną ... piękne wzory tylko nie da się tego policzyć !!! ;)
- 4 -
Reprezentacja w przestrzeni
ℝ
N
W a ,b =
〈
x n ,
a , b
n
〉
=
∑
n=0
N
x n
a , b
n
Ważne definicje:
–
położenie i rozciągłość w czasie
t =
∫
ℝ
t∣
a ,b
t ∣dt
t
2
=
∫
ℝ
t−t
2
∣
a , b
t∣dt
–
położenie i rozciągłość w częstotliwości
=
∫
ℝ
∣
a ,b
∣
d
2
=
∫
ℝ
−
2
∣
a ,b
∣
d
(rysunek z
t
i
)
Zasada nieoznaczoności Heisenberga:
t
≥
1
2
Przypadek dyskretny
Dwa rozwiązania:
1) nowe współrzędne skali
a
k
, l a
k
b gdzie l , k ∈ℤ
k , l
n=a
−
k / 2
n−l b
a
m
2) współrzędne skali z podziałem przez 2
a
k
=
a
0
k
, a
0
=
2 a
k
=
2
k
b
l
=
l b
0
a
0
k
, b
0
=
1 b
l
=
l2
k
- 5 -
k , l
n=
1
a
k
n−b
l
a
k
=
2
−
k /2
2
−
k
n−l
Zatem w przestrzeni
ℝ
N
W k , l=
〈
x n ,
k ,l
n
〉
=
∑
n=0
N −1
x n
k , l
n
i konsekwentnie transformata odwrotna:
x n=
∑
k
∑
l
W k , l
k ,l
n
Falka Morleta
Prototyp
n=
1
2
e
−
j
0
n
e
−
n
2
/
2
=
1
2
e
−
j
0
n n
2
/
2
=
e
−
2
2
−
0
2
zwykle
0
=
2
ln 2
=
5.336
Meksykański kapelusz
n=1−n
2
e
−
n
2
/
2
N=4;n=(-N:.1:N);
psi=(1-n.^2).*exp(-(n.^2)./2);
psi_1_0=sqrt(2)*(1-(n/2).^2).*exp(-((n/2).^2)./2);
- 6 -
plot(n,real(psi),';Re;',n,imag(psi),';Im;');
- 7 -