Zaawansowane algorytmy analizy częstotliwościowej sygnałów
Autorzy: Tomasz P. Zieliński
Pliki: rozdzial_20.zip
Przykład 20.1
W przykładzie 20.1 (plik przyklad_20_1.m) przedstawiono w języku MATLAB dekompozycję EVD macierzy funkcji autokorelacji sygnału. Jest w nim generowana suma K, od 1 do 4, nietłumionych lub tłumionych sygnałów sinusoidalnych, zaszumonych. Należy zmieniać liczbę składowych, wartości ich pulsacji i tłumień, poziom szumu i obserwować wynik dekompozycji.% Przykład 20.1
% Autor: Tomasz ZIELIŃSKI, tzielin@agh.edu.pl
clear all; close all;
%% Parametry sygnałów
N = 500; % liczba próbek sygnału
K = 2; % liczba składowych rzeczywistych
SNR = 160; % stosunek sygnał/szum [decybele], od 20 do 160
Om = [12.345 77.89 111.22 145.67 ]*(2*pi)/N, % pulsacje - numery prążków DFT
d = [0.001 0.01 0.10 0.20 ], % tłumienia
% WAŻNE: dla fPisarenko i fESPRIT wymagane d=0
%d = [ 0 0 0 0 ]; % dla metod fPisarenko i fESPRIT
A = [4.3 3.3 2.3 1.3 ]; % amplitudy
p = [-pi/11 -pi/7 -pi/5 -pi/3]; % fazy
%% Generacja sygnału testowego
n = 0:N-1; x = zeros(1,N);
for k=1:K
x = x + A(k)*cos(Om(k)*n+p(k)).*exp(-n*d(k));
end
y = awgn( x, SNR); % dodanie szumu
plot(y); title('Analizowany sygnał'); xlabel('n'); grid; pause
%% Dekompozycja EVD
L = 2*K+20; % wymiar macierzy autokorelacji
ryy=xcorr(y,'unbiased'); % do wyboru: ‘unbiased', ‘biased', 'coeff', 'none'
ryy = ryy(N:end); % z powodu symetrii, bierzemy tylko połowę
Ryy = toeplitz( ryy(1:L+1) ); % zbudowanie macierzy Rxx
[U, D ] = eig( Ryy ); % EVD( Rxx )
sigma = diag( D ); % wartości własne
[sigma, indx] = sort( sigma, 'descend' );
stem( sigma ), title('Warości własne'); grid, pause
subplot(411); plot(U(:,indx(1))), title('Wektor własny 1'); grid,
subplot(412); plot(U(:,indx(2))), title('Wektor własny 2'); grid,
subplot(413); plot(U(:,indx(3))), title('Wektor własny 3'); grid,
subplot(414); plot(U(:,indx(4))), title('Wektor własny 4'); grid,
pause
Przykład 20.2
W przykładzie 20.2 (plik fHilbert.m)zastosowano funkcję środowiska MATLAB hilbert() i pokazano, jak można wykorzystać transformację Hilberta do estymacji pulsacji i tłumienia pojedynczej składowej sinusoidalnej (20.1). Obliczane są w nim chwilowe wartości tych wielkości dla każdej próbki sygnału, a następnie są one uśredniane. Ponieważ na początku i na końcu okna obserwacji występują wyraźne oscylacje w obu estymatach (pulsacji i tłumienia), dlatego odrzucane są obliczone N/4 próbki początkowe i końcowe obu estymat. W ostatniej linii programu użyta jest zmienna temp, która służy do estymacji wartości tłumienia metodą regresji liniowej (aproksymacja linią prostą punktów pomiarowych). Sygnał, który będzie analizowany, można wygenerować w następujący sposób:N=1000; fpr=1000; dt=1/fpr; t=dt*(0:N-1); dA = 0.5; fA=2; A = 10*(1+dA*sin(2*pi*fA*t)); df=200; fn = 250; fm=4; fi=2*pi*(fn*t-df/(2*pi*fm)*cos(2*pi*fm*t)); x= A.*cos(fi);Jest to sygnał z równoczesną modulacją amplitudy i częstotliwości. Na początku przyjmujemy dA = 0, df=0, a potem dA = 0.5, df=200. W drugim przypadku usuwamy znak % przed instrukcjami plot().
function [Ome, de]=fHilbert(x,W) % Przykład 20.2 % x = A*cos(Om*n+p).*exp(-d*n) % W - parametr wyboru okna: 0-Hanning, 1-Blackman % Autor: Tomasz ZIELIŃSKI, tzielin@agh.edu.pl N = length(x); % długość sygnału if(W==1) win = hanning(N)'; end % wybór okna if(W==2) win = blackman(N)'; end % Xh = hilbert(win.*x)./win; % transformacja Hilberta w dziedzinie częstotliwości Om = diff( unwrap( angle(Xh) ) ); % plot(Om), pause % pulsacja chwilowa - zmiany A = abs(Xh); % plot(A), pause % amplituda chwilowa - zmiany d = -diff( log(A) ); % plot(d), pause % tłumienie chwilowe - zmiany ind = round( N/2-N/4 : N/2+N/4 ); % obcięcie oscylacji na brzegach Ome = mean(Om(ind)); % uśredniona wartość pulsacji de = mean(d(ind)); % uśredniona wartość tłumienia % temp = polyfit(ind,log(abs(Xh(ind))),1); de=-temp(1); % regresja liniowa
Przykład 20.3
W przykładzie 20.3 (plik fProny.m), napisanym w języku MATLAB, zaimplementowano opisaną powyżej metodę, opartą na równaniach (20.41), (20.44), (20.49c). Zastosowana wewnętrzna funkcja MATLABa toeplitz() jest w nim wykorzystana do zbudowania macierzy X, takiej jak w (20.48) – jako parametry wejściowe otrzymuje ona pierwszą kolumnę (x(2:N-1)) oraz pierwszy wiersz (x(2:-1:1)) macierzy X. Zwróćmy uwagę, że MATLAB indeksuje elementy wektorów i macierzy od 1, a nie od 0. Do funkcji możemy przekazać do analizy sygnał, generowany w przykładzie 20.1.function [Ome, de]=fProny(x,K)
% Przykład 20.3
% Metoda Prony'ego
% x - analizowany sygnał - suma tłumionych sinusoid rzeczywistych
% K - założona liczba składowych
% Autor: Tomasz ZIELIŃSKI, tzielin@agh.edu.pl
% DLA JEDNEJ SKŁADOWEJ)
if(K==1)
N = length(x); X = toeplitz(x(2:N-1), x(2:-1:1)); x = x(3:N).'; % (20.48)
a = X\x; % lub a = pinv(X)*x lub inv(X'*X)*X'*x;; % (20.49c)
p = log( roots( ([1; -a]) )); % (20.43)
Ome = abs(imag(p(1))); de = -real(p(1)); % (20.44)
end
% DLA KILKU SKŁADOWYCH
if(K>1)
N = length(x);
P = 2*K; % rząd predykcji
FirstCol = x(P:N-1); FirstRow = x(P:-1:1); %
X = toeplitz( FirstCol, FirstRow ); x = x(P+1:N).'; % (X.48)(X.57)
a = X\x; % lub inv(X'*X)*X'*y; lub a = pinv(X)*y; % (X.49c)(X.58)
a1 = [1; -a]; % rozwiązanie Prony’ego
p = roots( a1 ); %(X.43)(X.59)
p = log(p) ; % (X.44)
Om = imag(p); [Om indx] = sort( Om, 'descend' ); Ome = Om(1:K); % X.44)
de = -real(p(indx(1:K)));
end
Przykład 20.4
W przykładzie 20.4 (plik fSTMB.m), zawierającym implementację algorytmu Steiglitza–McBride’a w języku MATLAB, algorytm ten jest kodowany z uwzględnieniem wzorów (20.56)–(20.63) dotyczących sygnału, który może się składać z większej liczby tłumionych sinusoid (K >= 1).Wprzypadku przyjęcia zerowej liczby iteracji, wyliczane jest jedynie rozwiązanie równania (20.57), czyli wynik metody Prony’ego. Ta część programu stanowi uogólnienie kodu zaprezentowanego w przykładzie 20.3. Do funkcji możemy przekazać do analizy sygnał generowany w przykładzie 20.1.function [Ome, de]=fSTMB(x,K,NI) % Przykład 20.4 % x - analizowany sygnał % K - liczba tłumionych sinusoid A*cos(Om*n+p).*exp(-d*n) w równaniu (20.3) % NI - liczba iteracji, dla NI=0 otrzymujemy metodę Prony'ego % Ome - estymata pulsacji (rd) % de - estymata tłumienia % Autor: Tomasz ZIELIŃSKI, tzielin@agh.edu.pl N = length(x); % liczba próbek sygnału P = 2*K; % rząd predykcji FirstCol = x(P:N-1); FirstRow = x(P:-1:1); % pierwsza kolumna i rząd macierzy X (20.57) X = toeplitz( FirstCol, FirstRow ); y = x(P+1:N).'; % macierz X w (20.57), (20.48) a = X\y; % lub inv(X'*X)*X'*y; lub a = pinv(X)*y; % rozwiązanie (20.58), (20.49c) a1 = [1; -a]; % inicjalizacja rozwiązaniem metodą Prony'ego %a1 = [1; zeros(P,1) ]; % prostsza inicjalizacja zerami delta = zeros(1,N); delta(1)=1; % impuls jednostkowy while ( NI ) v = filter(1, a1, x); % (20.51) u = filter(1, a1, delta); % (20.52) V = toeplitz( v, [ v(1) zeros(1,P) ] ); % (20.54) (20.63) U = toeplitz( u, [ u(1) zeros(1,P-1) ] ); % (20.54) (20.63) VU = [ V(1:N,:) U(1:N,:) ]; % (20.54) (20.63) A = VU(:,2:end); d = VU(:,1); % (20.54) (20.63) c = A\d; % c = inv(A'*A)*A'*d; % c = pinv(A)*d; % (20.55) a1 = [ 1; -c(1:P) ]; % iteracyjna aktualizacja NI = NI-1; end p = roots( a1 ); % (20.43)(20.59) p = log(p); % (20.44) Om = imag(p); [Om indx] = sort( Om, 'descend' ); Ome = Om(1:K); % (20.44) de = -real(p(indx(1:K))); % (20.44)
Przykład 20.5
W przykładzie 20.5 (plik flpCOR.m) przedstawiono kod w języku MATLAB, pokazujący jak metody korelacyjne implementuje się programowo. Do funkcji możemy przekazać do analizy sygnał generowany w przykładzie 20.1.function [Ome, de] = flpCOR(x,K,nr) % Przykład 20.5 % x - analizowany sygnał % K - liczba tłumionych sinusoid A*cos(Om*n+p).*exp(-d*n) w równaniu (20.3) % nr - numer metody % Autor: Tomasz ZIELIŃSKI, tzielin@agh.edu.pl N = length(x); % liczba próbek sygnału L = 2*K; % rząd predykcji, podwojona liczba składowych rzeczywistych M = N-L; % maksymalna liczba równań r = xcorr(x,'unbiased'); % funkcja korelacji: ‘unbiased', ‘biased', 'coeff', 'none' r = r(N:end); % tylko przesunięcia nieujemne w funkcji korelacji if(nr == 1) % ŚREDNIO DOBRE, GRORSZE NIŻ NR=2 R1 = toeplitz( r(1:L) ); % macierz LxL r1 = r(2:L+1).'; % pause % macierz Lx1 a = R1\r1; end if(nr == 2) % ZAWSZE BARDZO DOBRE FirstCol = r(L:L+L-1); FirstRow = r(L:-1:1); R2 = toeplitz( FirstCol, FirstRow ); % macierz LxL r2 = r(L+1:L+L).'; % pause % macierz Lx1 a = R2\r2; end if(nr == 3) % DOBRE TYLKO BEZ TŁUMIEŃ FirstCol = r(L:M+L-1); FirstRow = r(L:-1:1); R3 = toeplitz( FirstCol, FirstRow ); % macierz MxL r3 = r(L+1:L+M).'; % pause % macierz Mx1 a = R3\r3; end a1 = [1 -a.']; p = roots( a1 ); p = log(p) ; % (20.43)(20.44)(20.59) Om = imag(p); [Om indx] = sort( Om, 'descend' ); Ome = Om(1:K); % (20.44) de = -real(p(indx(1:K)));
Przykład 20.6
W przykładzie 20.6 (plik fKT.m) jest przedstawiony kod w języku MATLAB funkcji fKT(), w której zaimplementowano algorytm Kumaresana–Tuftsa. Jest to zmodyfikowana wersja procedury lpsvd.m, pochodzącej z biblioteki matNMR [4]. W przeciwieństwie do oryginału działa ona na sygnałach o wartościach rzeczywistych, a nie na zespolonych. Do funkcji możemy przekazać do analizy sygnał, generowany w przykładzie 20.1.function [Ome, de]=fKT(x,K) % Przykład 20.6 % x - sygnał analizowany % K - założona liczba składowych (tłumione sinusoidy o~wartości rzeczywistej) % Autor: Tomasz ZIELIŃSKI, tzielin@agh.edu.pl (na podstawie lpsvd.m [4]) M = 2*K; % liczba składowych zespolonych N = length(x); % liczba próbek L = floor(N*3/4); % przyjęty rząd liniowej predykcji L = 3/4*N Y = hankel( x(2:N-L+1), x(N-L+1:N) ); % macierz próbek, związana z predykcją "do tyłu" y = x(1:N-L)'; % wektor próbek, związana z predykcją "do tyłu" [U,S,V] = svd(Y,0); % dekompozycja względem wartości szczególnych S = diag(S); szum = mean(S(M+1:min([N-L,L]))); % przesunięcie - średnia energia szumu b = V(:,1:M) * (diag(1./(S(1:M)-szum)) * (U(:,1:M)'*y)); % współczynniki predykcji "do tyłu" z = roots([-b(length(b):-1:1);1]); p=log(z); % pierwiastki wielomianu (20.73), ln() p = p( find(real(p)<0) ); % wewnątrz okręgu jednostkowego Om = imag(p); [Om indx] = sort( Om, 'descend' ); Ome = Om(1:K); % pulsacja kątowa de = -real(p(indx(1:K))); % współczynnik tłumienia
Przykład 20.7
W przykładzie 20.7 (plik fMatPen.m) zaprezentowano funkcję fMatPen() [36] napisaną w języku MATLAB. Kod w języku Fortran można znaleźć w [29]. Do funkcji możemy przekazać do analizy sygnał, generowany w przykładzie 20.1.function [Ome, de] = fMatPen(x,K) % Przykład 20.7 % x - analizowany sygnał - suma tłumionych sinusoid rzeczywistych % K - założona liczba składowych % Autor: Tomasz ZIELIŃSKI, tzielin@agh.edu.pl M = 2*K; % liczba składowych zespolonych N = length(x); % liczba próbek sygnału L = floor(N/3); % przyjęty rząd liniowej predykcji L = N/3 X = hankel(x(1:N-L),x(N-L:N)); % X0=X(:,1:L), X1=X(:,2:L+1) [U,S,V] = svd(X(:,2:L+1), 0); S = diag(S); % SVD( X1 ) p = log( eig( diag(1./S(1:M)) * ((U(:,1:M)'*X(:,1:L))*V(:,1:M)) ) ); % ln(EVD( Z )) Om = imag(p); [Om indx] = sort( Om, 'descend' ); Ome = Om(1:K); % pulsacja de = real(p(indx(1:K))); % tłumienie
Przykład 20.8
W przykładzie 20.8 (plik fPisarenko.m) jest zaimplementowana w języku MATLAB jedynie metoda Pisarenki, ponieważ bezpośrednio można w niej wyliczyć poszukiwane wartości parametrów sygnału. W innych metodach znajduje się je metodą poszukiwania maksimów odpowiedniej funkcji kosztu, co jest bardziej kłopotliwe. Nie jest to duża strata, gdyż brakujące algorytmy mają gorszą skuteczność, niż opisane metody Kumaresana–Tuftsa, TLS i Matrix Pencil. Do funkcji możemy przekazać do analizy sygnał generowany w przykładzie 20.1, z d = [0, 0, 0, 0].function [Ome]=fPisarenko(y,K) % Przykład 20.8 % y - analizowana suma tłumionych sinusoid o wartościach rzeczywistych w szumie % K - założona liczba składowych % Autor: Tomasz ZIELIŃSKI, tzielin@agh.edu.pl N = length(y); % liczba próbek sygnału L = 2*K; % rząd predykcji, liczba składowych zespolonych ryy=xcorr(y,'unbiased'); % do wyboru: 'unbiased', 'biased', 'coeff', 'none' ryy = ryy(N:end); % z powodu symetrii, bierzemy tylko połowę wyniku Ryy = toeplitz( ryy(1:L+1) ); % zbudowanie macierzy Ryy [V,D] = eig(Ryy); % dekompozycja EVD tej macierzy [Dmin indx] = min(diag(D)); % znalezienie najmniejszej wartości własnej a = V(:,indx).'; % znalezienie skojarzonego z nią wektora własnego z = roots(a); Om = imag(log(z)); % znalezienie pierwiastków wielomianu Om = sort( Om, 'descend' ); Ome = Om(1:K); % pulsacja (odrzucenie wartości ujemnych)
Przykład 20.9
W przykładzie 20.9 (plik fESPRIT.m) jest przedstawiony kod w języku MATLAB algorytmu ESPRIT. Do funkcji możemy przekazać do analizy sygnał generowany w przykładzie 20.1, z d = [0, 0, 0, 0].function [Ome]=fESPRIT(x,K) % Przykład 20.9 % x - analizowany sygnał - suma tłumionych sinusoid o~wartościach rzeczywistych % K - założona liczba składowych % Autor: Tomasz ZIELIŃSKI, tzielin@agh.edu.pl N = length(x); % licza próbek sygnału L = 2*K; % rząd predykcji, liczba składowych zespolonych L1 = L+1; r=xcorr(x,'unbiased'); % do wyboru: 'unbiased', 'biased', 'coeff', 'none' r = r(N:end); % funkcja autokorelacji ryy R0 = toeplitz( r(1:L+1) ); % macierz autokorelacji Ry0y0 R1 = [ R0(:,2:L1) r(L1+1:-1:2)']; % macierz autokorelacji Ry0y1 [V,D] = eig(R0); % EVD( Ry0y0 ) sigma = min(diag(D)); % minimalna wartość własna I = eye(L1); J = [zeros(1,L1); I(1:L1-1,1:L1) ]; % macierze I, J C0 = R0 - sigma*I; C1 = R1 - sigma*J; % macierze C0= Cy0y0. C1=Cy0y1 S = eig(C0, C1); % EVD( C0, C1) - Matrix-Pencil S = S( find( abs(S) > 0.5 ) ); % znajdź wartości zespolone o module > 0.5 Om = atan2( imag(S), real(S) ); % kąt liczby zespolonej Om = sort( Om, 'descend' ); Ome = Om(1:K); % pulsacje (odrzucenie wartości ujemnych)