Gdańsk, 01.10.2005

METODY IDENTYFIKACJI - PROJEKT

Część I: Zastosowanie metod identyfikacji do analizy i syntezy sygnałów - opis tematów

Część II: Przetwarzanie sygnałów - programy narzędziowe

___________________________________________________________________

Część I - opis tematów

___________________________________________________________________

Cel projektu

Projekt ma na celu zapoznanie studentów z praktycznymi aspektami identyfikacji procesów. Od uczestników zajęć wymaga się znajomości metod i algorytmów identyfikacji omawianych w ramach wykładu „Identyfikacja procesów” prowadzonego na III roku studiów na kierunku Automatyka i Robotyka. Zadania projektowe dotyczące predykcyjnego kodowania sygnałów, identyfikacji strukturalnej modeli z czasem dyskretnym oraz estymacji parametrycznej ciągłych modeli obiektów stanowią ilustrację tzw. parametrycznego podejścia do analizy procesów, wykorzystywanego powszechnie w nowoczesnych systemach automatyki i telekomunikacji.

W ramach projektu praktycznego (predykcyjne kodowanie sygnałów) badania przeprowadzane są na rzeczywistych sygnałach (nagraniach monofonicznych), co pozwala ocenić uzyskane wyniki nie tylko obiektywnie (tj. w oparciu o wartości odpowiednich wskaźników jakości), ale i subiektywnie (tzn. odsłuchowo). Realizując z kolei projekty teoretyczne (identyfikacja strukturalna modeli z czasem dyskretnym i estymacja parametryczna ciągłych modeli obiektów) uczestnicy zajęć mają okazję zweryfikować w numerycznym środowisku symulacyjnym (MATLAB) działanie algorytmów omawianych na wykładzie.

Warunkiem zaliczenia przedmiotu jest wykonanie co najmniej dwóch zadań.

Opis stanowiska - sprzęt i narzędzia

Do dyspozycji uczestników zajęć przewidziano stanowiska wyposażone w komputery osobiste z zainstalowanym środowiskiem symulacyjnym MATLAB. Dodatkowo uczestnicy zajęć mogą wykorzystać programy narzędziowe (view.exe, wave2snd.exe, snd2wave.exe) opisane w części II. Przy rozliczaniu projektów dostępny będzie również odpowiedni sprzęt nagłaśniający oraz słuchawki (odsłuchowa weryfikacja jakości nagrań).

Dopuszcza się realizację zadań projektowych z wykorzystaniem innych narzędzi programistycznych (np. kompilatorów języka Pascal lub języka C/C+ + ), jeśli tylko uczestnik zajęć posiada (legalnie) takie oprogramowanie.

Postać sprawozdania

Dla każdego zrealizowanego zadania dwuosobowa grupa projektowa wykonuje osobne sprawozdanie zawierające:

a) Krótki opis podejścia wraz z informacją o parametrach projektowych i własnościach programu.

b) Charakterystykę otrzymanych wyników.

c) Wydruk programu źródłowego.

d) Wnioski końcowe.

Napisane programy powinny zawierać wygodny dla użytkownika interfejs.

Projekt 1 - predykcyjne kodowanie sygnałów

W pliku proj1_nr.wav dany jest ok. 10 sekundowy fragment zapisu sygnału mowy (każda próbka reprezentowana przez 16-bitowe słowo, częstotliwość próbkowania 11025 Hz). Napisać program (nadajnika) kodujący zarejestrowany sygnał przy użyciu odpowiedniego filtru 10 rzędu oraz program (odbiornika) rekonstruujący sygnał na podstawie przesłanych z nadajnika parametrów filtru.

Procedura nadajnika

a) Podzielić zarejestrowany sygnał na segmenty obejmujące N próbek każdy (N = 256), przy czym sąsiednie segmenty nakładają się (tj. zawierają r (r = 10) tych samych próbek na krańcach)

b) Sygnał w każdym segmencie zamodelować równaniem autoregresyjnym (AR) rzędu r (r = 10)

, k = 1 , ... , N

gdzie y(k) oznacza wartość próbki sygnału, a_i (i = 1 , ... , r) są współczynnikami zastosowanego filtru 1/A(z ^-1) = 1/( 1 + a₁ z ^-1 + ... + a_r z ^-r ), zaś e(k) reprezentuje błąd resztowy. Z uwagi na założenie o stabilności filtru 1/A(z ^-1) współczynniki równania AR(r) zidentyfikować przy użyciu algorytmu Levinsona-Durbina (L-D), a nie klasycznego algorytmu najmniejszych kwadratów. Procedurę L-D należy napisać i zaimplementować samodzielnie.

Przed uruchomieniem algorytmu identyfikacji dane { y(1) , ... , y(N)} w obrębie segmentu trzeba odpowiednio „spłaszczyć” na krańcach. W tym celu wartości próbek należy pomnożyć przez pewne wagi w(k). Odpowiednie okno formujące można przyjąć jako

, k = 1 , ... , N

Uformowany w ten sposób segment należy uzupełnić po obydwu stronach ciągami r zer, uzyskując dane do procedury identyfikacji: {0 , ... , 0 , w(1) y(1) , ... , w(N) y(N) , 0 , ... , 0}.

c) Wyznaczyć ciąg {e(k)} błędów resztowych: e(k) = y(k) + a₁ y(k - 1) + ... + a_r y(k - r) dla k = 1 , ... , N. W obrębie segmentu znaleźć największy (co do wartości bezwzględnej) błąd resztowy (e_max), a następnie przeprowadzić równomierną m-bitową kwantyzację wszystkich błędów w przedziale wartości <- e_max , e_max 〉.

d) Współczynniki a_i (i = 1 , ... , r) filtru, wartość e_max oraz ciąg {
} skwantowanych błędów resztowych przesłać do odbiornika. Procedurę powtórzyć dla kolejnych segmentów danych.

Procedura odbiornika

a) Na podstawie otrzymanych współczynników a_i (i = 1 , ... , r) filtru, wartości e_max oraz ciągu {
} skwantowanych błędów resztowych zrekonstruować sygnał w danym segmencie

, k = 1 , ... , N

Jako warunki początkowe { y(0) , ... , y(1 - r)} przyjąć r próbek poprzedniego segmentu.

b) Odtworzony segment danych dopisać do pliku zawierającego zrekonstruowane próbki nadawanego sygnału. Procedurę powtórzyć dla kolejnych danych otrzymywanych z nadajnika.

Napisany program nadajnika-odbiornika należy uruchomić dla różnych poziomów kwantyzacji błędów resztowych, tj. dla m = 2 (kwantyzacja dwubitowa - 4-poziomy), dla m = 3 (kwantyzacja trzybitowa - 8-poziomów) i dla m = 4 (kwantyzacja czterobitowa - 16-poziomów). Plik oryginalny oraz pliki zrekonstruowane należy porównać używając edytora view.exe, zaś jakość rekonstrukcji należy ocenić odsłuchowo.

Literatura: Spanias A. (1994): „Speech coding: a tutorial review”. Proc. of the IEEE, vol.82, no.10, str.1541-1562.

Projekt 2 - identyfikacja strukturalna modeli z czasem dyskretnym

W pliku proj2_nr.bin dany jest ciąg N (N = 4000) próbek { y(1) , ... , y(N)} uzyskanych z równania autoregresyjnego AR(r) o nieznanym rzędzie r (3 ≤ r ≤ 12)

, k = 1 , ... , N

gdzie e(k) jest gaussowskim szumem białym o pewnej wariancji. Każda próbka jest zapisana jako liczba rzeczywista 8-bajtowa (typ float64 w środowisku MATLAB lub typ double w języku C). Napisać program odczytujący dane z pliku oraz dokonujący strukturalnej i parametrycznej identyfikacji modelu AR(r). Po uzyskaniu rzędu r oraz parametrów a_i (i = 1 , ... , r) model należy zweryfikować.

Procedura identyfikacji strukturalnej i weryfikacji

a) Zapisać równanie autoregresyjne w postaci wektorowej (k = 1 , ... , N )

z wektorami regresji i parametrów zdefiniowanymi następująco

b) Stosując klasyczny algorytm najmniejszych kwadratów (LS) i przyjmując stałą zapominania λ = 1 wyznaczyć oceny parametrów modelu AR(r) dla rzędów r = 3 , ... , 12

, P(0) = diag (10⁵ , ... , 10⁵ )

,

oszacowując każdorazowo wariancję σ_e² błędu resztowego

,

c) Wykorzystując kryterium końcowego błędu predykcji (FPE) oraz informacyjne kryteria Akaike'go (AIC) ocenić właściwy rząd modelu. Rozważane funkcje kryterialne

,
,

wykreślić w funkcji parametru r wskazując na wykresach tę wartość r^*, dla której rozpatrywane wskaźniki osiągają minimum. Zapisać współczynniki a_i (i = 1 , ... , r) zidentyfikowanych modeli rzędów r = r^*, r = r^* - 1 i r = r^* - 2.

d) Dokonać weryfikacji modelu AR(r^*) z parametrami
wykreślając unormowaną funkcję autokorelacji γ (m)

,

i sprawdzając, czy spełnione są warunki

, m = 1 , ... , N - 1

Adekwatność modelu zbadać także wykreślając periodogram kumulacyjny g(m)

,

i sprawdzając, czy zachodzą nierówności

, m = 1 , ... , [N / 2]

W obliczeniach zastosować η_e = 1.36, co odpowiada zalecanemu poziomowi ufności 5%.

e) Obliczyć ze wzoru definicyjnego widmową gęstość mocy S(ω) procesu y(k)

,

a następnie, korzystając z wyników identyfikacji, wyznaczyć parametryczną ocenę widma mocy

Uzyskane funkcje S(ω) i
przedstawić na wspólnym wykresie.

f) Korzystając z otrzymanego modelu AR(r^*) wygenerować nowy ciąg N (N = 4000) próbek szeregu czasowego { y(k)}. Zastosować generator gaussowskiego szumu białego przyjmując wariancję na poziomie σ_e² (w środowisku MATLAB wykorzystać funkcję randn( 1 , N ) ).

W programie powinny się znaleźć procedury umożliwiające wykreślanie na ekranie funkcji kryterialnych FPE(r) i AIC(r), funkcji autokorelacji γ (m), periodogramu g(m) oraz widma mocy S(ω) i oceny widma
.

Literatura: Box G.E.P., Jenkins G.M. (1983): „Analiza szeregów czasowych. Prognozowanie i sterowanie”. PWN, Warszawa.

Projekt 3 - estymacja parametryczna modeli z czasem ciągłym

W plikach proj3u_nr.bin i proj3y_nr.bin dane są N-elementowe (N = 3600) ciągi próbek u(k) = u(kT) i y(k) = y(kT) sygnału u(t) podawanego na wejście obiektu i sygnału y(t) obserwowanego na wyjściu obiektu (okres próbkowania T = 0.05 s, czas symulacji t_max = 180 s). Zakłada się, że dynamikę obiektu opisuje liniowe równanie różniczkowe zwyczajne rzędu r (r = 3)

zaś pomiar y(k) sygnału y(t) jest zakłócony gaussowskim szumem białym v(k) o pewnej wariancji: y(k) = y(t) _| t = kT + v(k). Korzystając z pośredniej i bezpośredniej metody identyfikacji modeli z czasem ciągłym oraz stosując odpowiednie algorytmy estymacji parametrycznej wyznaczyć oceny parametrów α_i i β _i (i = 1 , ... , r ) modelu.

Procedura estymacji parametrycznej - metoda pośrednia

a) Przyjąć pomocniczy model autoregresji z zewnętrznym wejściem (ARX) rozważanego obiektu ciągłego, tj. y(k) + a₁ y(k - 1) + ... + a_r y(k - r) = b₁ u(k - 1) + ... + b_r u(k - r) + e(k). Równanie ARX zapisać w postaci wektorowej (k = 1 , ... , N )

z wektorami regresji i parametrów zdefiniowanymi następująco

b) Stosując klasyczny algorytm najmniejszych kwadratów (LS) i przyjmując stałą zapominania λ = 1 wyznaczyć oceny parametrów a_i i b_i (i = 1 , ... , r ) modelu ARX

, P(0) = diag (10⁵ , ... , 10⁵ )

,

c) Wykorzystując oceny LS parametrów a_i i b_i dokonać powtórnej identyfikacji modelu asymptotycznie nieobciążoną metodą zmiennej instrumentalnej (IV) ze stałą λ = 1

, P(0) = diag (10⁵ , ... , 10⁵ )

,

Nieskorelowany z zakłóceniem wektor zmiennych pomocniczych ξ(k) można przyjąć jako

lub

lub

przy czym w ostatnim przypadku ocenę niezakłóconego wyjścia obiektu dostaje się stosując w każdym kroku identyfikacji filtr

, k = 1 , ... , N

z ocenami
i
parametrów uzyskiwanymi każdorazowo z poprzedniego kroku algorytmu IV.

d) Dokonać „rzutowania” zidentyfikowanego modelu dyskretnego z płaszczyzny operatora przesunięcia z na płaszczyznę operatora różniczkowania s (Laplace'a). Rekonstruując model ciągły przyjąć założenie o niezmienniczości odpowiedzi skokowej (ZOH). Dopuszcza się wykorzystanie funkcji bibliotecznej pakietu MATLAB: Gs = d2c(Gz , 'zoh'), gdzie Gz oznacza znaną (zidentyfikowaną) z-transmitancję, zaś Gs jest poszukiwaną s-transmitancją.

Procedura estymacji parametrycznej - metoda bezpośrednia

a) Podane równanie różniczkowe scałkować obustronnie stosując operator całkowania ze skończonym horyzontem (L) obserwacji: J ^r(s) = [ ( 1 - e ^-sLT ) / s] ^r. Stosując podstawienie Tustina: s ^-1 = 0.5 T ( 1 + z ^-1 ) / ( 1 - z ^-1 ), dokonać dyskretyzacji opisu ciągłego otrzymując model w postaci wektorowej (k = 1 , ... , N )

z wektorami regresji i parametrów zdefiniowanymi następująco

Zastosowane filtry formujące typu FIR (r - krotne całki kolejnych pochodnych j = 0 , ... , r ) wyznaczyć ze wzoru ( L ∈ <20 , 40〉 )

b) Metodą najmniejszych kwadratów (λ = 1) znaleźć oceny parametrów α_i i β _i modelu ciągłego

, P(0) = diag (10⁵ , ... , 10⁵ )

,

c) Wykorzystując oceny LS parametrów α_i i β _i dokonać powtórnej identyfikacji modelu asymptotycznie nieobciążoną metodą zmiennej instrumentalnej ze stałą λ = 1

, P(0) = diag (10⁵ , ... , 10⁵ )

,

Nieskorelowany z zakłóceniem wektor zmiennych pomocniczych ξ(k) można przyjąć jako

przy czym ocenę niezakłóconego wyjścia obiektu dostaje się stosując w każdym kroku identyfikacji filtr

, k = 1 , ... , N

z ocenami
i
parametrów uzyskiwanymi każdorazowo z poprzedniego kroku algorytmu IV.

d) Wykreślić unormowany moduł transmitancji widmowej zidentyfikowanego obiektu oraz rodzinę krzywych reprezentujących unormowany moduł transmitancji widmowej filtru całkującego dla różnych wartości horyzontu całkowania (L = 20, 25, 30, 35 i 40)

,

Na podstawie wykresu dobrać taką wartość horyzontu L, aby pasmo przenoszenia obiektu i pasmo przenoszenia filtru całkującego pokrywały się dla ω ∈ < 0 , 2π / (LT ) 〉. Dla tak określonej wartości horyzontu L powtórzyć identyfikację modelu metodami LS i IV.

Porównać wyniki identyfikacji uzyskane metodą pośrednią i metodą bezpośrednią. W napisanym programie umieścić procedury umożliwiające wykreślanie na ekranie ciągów
i
oraz
i
ocen parametrów generowanych metodami identyfikacji LS i IV.

Literatura: Sagara S., Zhao Z.Y. (1990): „Numerical integration approach to on-line identification of continuous-time systems. Automatica, vol.26, no.1, str.63-74.

___________________________________________________________________

Część II - programy narzędziowe

___________________________________________________________________

Opisane programy powinny być uruchamiane w systemie MS-DOS. Poprawność działania w środowisku Windows nie jest gwarantowana.

1. Program konwersji plików wave2snd.exe (autor: Krzysztof Cisowski)

Pliki zawierające nagrania muzyczne zapisywane są w tzw. formacie WAVE. Do obróbki danych wygodniej jest posługiwać się plikami zapisanymi w formacie SOUND. Pliki takie pozbawione są nagłówka formatu WAVE i zawierają wyłącznie wartości próbek sygnałów.

Uruchamiając program wave2snd.exe konwersji z formatu WAVE do SOUND należy podać:

- nazwę pliku źródłowego WAVE (Source file name), np. mowa.wav ,

- nazwę pliku wynikowego SOUND (Dest file name), np. dane.snd .

2. Program konwersji plików snd2wave.exe (autor: Krzysztof Cisowski)

Program dopisuje do pliku typu SOUND nagłówek formatu WAVE. Otrzymany plik wynikowy nadaje się do odtwarzania w systemie Windows.

Uruchamiając program snd2wave.exe konwersji z formatu SOUND do WAVE należy podać:

- nazwę pliku źródłowego SOUND (Source file name), np. probki.snd ,

- nazwę pliku wynikowego WAVE (Dest file name), np. slowa.wav .

- częstotliwość próbkowania sygnałów (Sampling rate), np. 11 kHz,

- liczbę kanałów (Channels), tj. 1 - dla nagrania mono, 2 - dla nagrania stereo,

- liczbę bitów przypadających na jedną próbkę (Bits per sample), np. 16.

3. Program edycji próbek sygnałów view.exe (autor: Krzysztof Cisowski)

Program view.exe jest edytorem próbek sygnałów zapisanych w postaci ciągu liczb 16-bitowych ze znakiem. Sygnały mogą być zapisane w plikach dyskowych zarówno w tzw. postaci surowej RAW (pliki z rozszerzeniem `raw'), czyli bez jakiegokolwiek nagłówka, jak również w klasycznym formacie WAVE (pliki z rozszerzeniem `wav') stosowanym w systemie Windows. Ponieważ program umożliwia synchroniczne oglądanie dwóch plików dźwiękowych, wywołanie programu może wyglądać następująco:

D:\ MET_IDENTYFIK \ MOJ_KATALOG \ view.exe mowa1.wav mowa2.raw

Po uruchomieniu programu na ekranie pojawia się pasek menu oraz dwa okna sygnałowe oznaczone literami A i B. W wyróżnionych polach okien można odczytać istotne informacje zapisane we wskazanych plikach oraz następujące parametry wykresów sygnałów:

a) w górnej części okna odczytać można numer próbki przy lewej krawędzi ekranu, nazwę pliku, numer próbki przy prawej krawędzi ekranu oraz rozmiar pliku (FL),

b) po prawej stronie okna można odczytać

- pozycję kursora (t),

- wartość X(t) próbki wskazywanej przez kursor,

- mnożnik (*) skali wykresu sygnału,

- początek (B) zaznaczonego bloku sygnału (po wciśnięciu klawisza F5),

- koniec (E) zaznaczonego bloku sygnału (po wciśnięciu klawisza F6),

- wskaźniki (A lub B) aktywności okien.

Do obsługi okien sygnałowych stosowane są następujące klawisze:

↑ ↓ - zwiększanie/zmniejszanie mnożnika dwa razy,

Alt↑ Alt↓ - zwiększanie/zmniejszanie mnożnika do maksymalnej wartości,

PageUp - zwiększanie składowej stałej (256),

PageDn - zmniejszanie składowej stałej (-256),

Ctrl PageUp - zwiększanie składowej stałej (8),

Ctrl PageDn - zmniejszenie składowej stałej (-8),

0 - zerowanie składowej stałej,

← → - przesunięcie na wykresie (512 punktowym) o 1 piksel w lewo/prawo,

Alt← Alt→ - przesunięcie na wykresie (512 punktowym) o 10 pikseli w lewo/prawo,

Ctrl← Ctrl→ - przesunięcie wzdłuż całego wykresu pół okna w lewo/prawo,

Home - wykres początku sygnału,

End - wykres końca sygnału,

F3 - zwiększanie rozdzielczości o 10 (normalna rozdzielczość to 1 piksel/próbkę),

F4 - zmniejszenie rozdzielczości o 10,

Ctrl F3 - zwiększanie rozdzielczości o 100,

Ctrl F4 - zmniejszenie rozdzielczości o 100,

Alt F3 - zwiększanie rozdzielczości o 1000,

Alt F4 - zmniejszenie rozdzielczości o 1000,

S - synchronizacja/desynchronizacja kursorów w oknach A i B

(zgodność numerów próbek w oknach),

Ctrl S - pełna synchronizacja/desynchronizacja kursorów w oknach A i B

(kursor, mnożniki, składowe stałe),

F5 - ustawienie początku zaznaczanego bloku,

F6 - ustawienie końca zaznaczanego bloku,

Alt F5 - kasowanie zaznaczenia początku bloku,

Alt F6 - kasowanie zaznaczenia końca bloku,

A B - wskazywanie aktywnego okna,

Tab - przełączanie aktywnego okna ( z A na B lub odwrotnie),

F2 - wyświetlenie informacji o ilości wolnej pamięci w komputerze,

F10 - uaktywnienie paska menu programu.

Wybierając menu programu (klawisz F10) uaktywnia się pole Input/Output złożone z kilku opcji (rezygnacja z wykonania opcji następuje po naciśnięciu klawisza Escape):

Open file         - otwieranie pliku w aktywnym oknie, przy czym nazwę pliku można podać jawnie (podopcja Enter file name) lub wybrać nazwę pliku z listy plików znajdujących się w aktualnym katalogu (podopcja Choose file name); żądany wariant wskazuje się naciskając klawisz Tab,

RAW-WAVE  - zmiana pliku z formatu RAW na WAVE (wybór pliku odbywa się analogicznie jak w opcji Open file),

WAVE-RAW  - zmiana pliku z formatu WAVE na RAW (wybór pliku odbywa się analogicznie jak w opcji Open file),

QUIT               - zakończenie programu.

256

segment ( l + 1)

segment ( l)

segment ( l - 1)

10

10

256

256

---