USG/Doppler - Historia wersji

Tsteifer: /* Metoda dopplerowska impulsowa (Pulsed Wave) */

2016-07-29T08:30:12Z

Metoda dopplerowska impulsowa (Pulsed Wave)

Mlew o 07:41, 5 lip 2016

2016-07-05T07:41:26Z

Mlew o 11:58, 3 lip 2016

2016-07-03T11:58:05Z

Mlew o 11:57, 3 lip 2016

2016-07-03T11:57:03Z

Mlew o 11:35, 3 lip 2016

2016-07-03T11:35:47Z

Mlew o 11:32, 3 lip 2016

2016-07-03T11:32:38Z

Mlew o 11:27, 3 lip 2016

2016-07-03T11:27:11Z

Mlew o 11:25, 3 lip 2016

2016-07-03T11:25:34Z

Mlew o 11:24, 3 lip 2016

2016-07-03T11:24:30Z

Mlew o 11:22, 3 lip 2016

2016-07-03T11:22:24Z

@@ Linia 3: / Linia 3: @@
 <source lang = python>
-f0 = 4e6 # Częstotliwość nadawcza [Hz]
+f0 = 5.5e6 # Częstotliwość nadawcza [Hz]
-z_step = 1.8970189701897084e-05 # Odległość między kolejnymi punktami w głębokości na siatce rekonstruowanego obrazu [m]
+z_step = 3.079291762894534e-05 # Odległość między kolejnymi punktami w głębokości na siatce rekonstruowanego obrazu [m]
 PRF = 1000  # Pulse Repetition Frequency [Hz]

← poprzednia wersja		Wersja z 07:41, 5 lip 2016
Linia 1:		Linia 1:
−	=Metoda dopplerowska=	+	=Metoda dopplerowska impulsowa (Pulsed Wave)=
	Do dyspozycji mamy zestaw zrekonstruowanych danych RF ([http://www.fuw.edu.pl/~jarekz/FUW-FID-2016/raw-data/usg3_doppler_circle_10lh.npy link]) z 64 kolejnych nadań pod stałym kątem (zerowym) falą płaską o parametrach:		Do dyspozycji mamy zestaw zrekonstruowanych danych RF ([http://www.fuw.edu.pl/~jarekz/FUW-FID-2016/raw-data/usg3_doppler_circle_10lh.npy link]) z 64 kolejnych nadań pod stałym kątem (zerowym) falą płaską o parametrach:

← poprzednia wersja		Wersja z 11:58, 3 lip 2016
Linia 146:		Linia 146:
	Porównać mapy prędkości przed i po filtracji medianowej (dla kilku różnych rozmiarów okna filtracji).		Porównać mapy prędkości przed i po filtracji medianowej (dla kilku różnych rozmiarów okna filtracji).

−	~~</hr>~~	+	----

@@ Linia 21: / Linia 21: @@
 Demodulację możemy przeprowadzić zarówno na danych surowych RF, jak i na danych po rekonstrukcji. W naszym wypadku zastosujemy to drugie rozwiązanie.
-Na początku stworzymy funkcję dokonującą demodulacji każdego z obrazów, tj. dla każdego próbki o współrzędnych <math>(x,y)</math>
+Na początku stworzymy funkcję dokonującą demodulacji każdego z obrazów, tj. dla każdej próbki o współrzędnych <math>(x,y)</math>
 <math>IQ(x,y)=RF(x,y)\cdot e^{-i\cdot2\pi f_{0}t} </math>
@@ Linia 27: / Linia 27: @@
 gdzie <math>IQ</math> - sygnał po demodulacji; <math>RF</math> - sygnał przed demodulacją; <math>t</math> - czas odpowiadający momentowi akwizycji próbki z danej głębokości; możemy przyjąć uproszczone założenie, że:
 <math>t=2y/c </math>; jak widać, pierwszy wymiar (szerokość) jest w naszej procedurze nieistotny.
-Po demodulacji sygnał należy przefiltrować dolnoprzepustowo w zakresie pasma podstawowego (tzw. baseband). W naszym wypadku wystarczający powinien być filtr o częstotliwości odcięcia równej 5.55MHz.
+Po demodulacji sygnał należy przefiltrować dolnoprzepustowo w zakresie pasma podstawowego (tzw. baseband). Dla naszych danych wystarczający powinien być filtr o częstotliwości odcięcia równej 5.55MHz.
 ===Zadanie===
@@ Linia 36: / Linia 36: @@
 ==Estymator autokorelacyjny==
 ===Rys teoretyczny===
-Dotychczas traktowaliśmy rekonstruowane obrazy jako sygnały dwuwymiarowe. Teraz, do zmierzenia średniej prędkości obiektów poruszających się w danym punkcie, konieczna będzie analiza wielu następujących po sobie obrazów, które zostały zebrane ze stałym interwałem czasowym (tzw. PRF - Pulse Repetition Frequency). W poniższych rozważaniach pomijamy wymiary przestrzenne - całe rozumowania przeprowadzone są dla sygnału <math>s(i)</math> na ustalonej głębokości i szerokości. Należy mieć na uwadze, że w praktyce ciężko mówić o estymacji prędkości w punkcie, ponieważ w naszej procedurze estymacji uwzględniamy efektywnie informacje z pewnego obszaru pomiarowego (energia danego piksela zawiera informację o obiektach z pewnego obszaru).
+Dotychczas traktowaliśmy rekonstruowane obrazy jako sygnały dwuwymiarowe. Teraz, do zmierzenia średniej prędkości obiektów poruszających się w danym punkcie, konieczna będzie analiza wielu następujących po sobie obrazów, które zostały zebrane ze stałym interwałem czasowym (tzw. PRF - Pulse Repetition Frequency). W poniższych rozważaniach pomijamy wymiary przestrzenne - całe rozumowania przeprowadzone są dla sygnału <math>s(i)</math> na ustalonej głębokości i szerokości. Należy mieć na uwadze, że w praktyce ciężko mówić o estymacji prędkości w punkcie, ponieważ w naszej procedurze estymacji uwzględniamy efektywnie informacje z pewnego obszaru pomiarowego (energia danego piksela zawiera informację z pewnej objętości pomiarowej, co wynika z kształtu wiązki ultradźwiękowej).
 Średnią prędkość będziemy mogli obliczyć korzystając ze średniego przesunięcia dopplerowskiego <math>\Delta f</math> w oparciu o szkolny wzór na [[Fizyka_III/Fale_dźwiękowe#Efekt_Dopplera|częstotliwość Dopplera]]:
@@ Linia 100: / Linia 100: @@
 ==Prezentacja Kolor==
 Po estymacji średnich prędkości w obszarach odpowiadających punktom na siatce użytej do rekonstrukcji obrazu, możemy nałożyć taką mapę na obraz B-mode w celu uzyskania obrazu "Kolor". Mając tablicę z danymi B-mode oraz tablicę prędkości, możemy otrzymać połączony obraz za pomocą poniższego skryptu:
 <source lang=python>
 # BMode - tablica zrekonstruowanych danych po przefiltrowaniu, obwiedni itp.
 # flow - tablica przepływów
@@ Linia 144: / Linia 145: @@
 ====Zadanie====
 Porównać mapy prędkości przed i po filtracji medianowej (dla kilku różnych rozmiarów okna filtracji).

@@ Linia 1: / Linia 1: @@
 =Metoda dopplerowska=
 Do dyspozycji mamy zestaw zrekonstruowanych danych RF ([http://www.fuw.edu.pl/~jarekz/FUW-FID-2016/raw-data/usg3_doppler_circle_10lh.npy link]) z 64 kolejnych nadań pod stałym kątem (zerowym) falą płaską o parametrach:
 <source lang = python>
 f0 = 4e6 # Częstotliwość nadawcza [Hz]
@@ Linia 26: / Linia 27: @@
 gdzie <math>IQ</math> - sygnał po demodulacji; <math>RF</math> - sygnał przed demodulacją; <math>t</math> - czas odpowiadający momentowi akwizycji próbki z danej głębokości; możemy przyjąć uproszczone założenie, że:
 <math>t=2y/c </math>; jak widać, pierwszy wymiar (szerokość) jest w naszej procedurze nieistotny.
-Po demodulacji sygnał należy przefiltrować dolnoprzepustowo w zakresie pasma podstawowego (tzw. baseband). W naszym wypadku wystarczający powinien być filtr dolnoprzepustowy o częstotliwości odcięcia równej 5.55MHz.
+Po demodulacji sygnał należy przefiltrować dolnoprzepustowo w zakresie pasma podstawowego (tzw. baseband). W naszym wypadku wystarczający powinien być filtr o częstotliwości odcięcia równej 5.55MHz.
-<b>Do wykonania:</b>
+===Zadanie===
 # Zbadać widmo amplitudowe sygnału przed i po demodulacji (po filtrowaniu). Jak zmieniło się widmo? Gdzie jest położona średnia widma? Czy rozkłady dla ujemnych i dodatnich częstotliwości są swoimi odbiciami?
 # Porównać widmo amplitudowe sygnału zdemodulowanego przed i po filtracji. Jakich składowych w częstości się pozbyliśmy (nie licząc szumu)?

@@ Linia 22: / Linia 22: @@
 Na początku stworzymy funkcję dokonującą demodulacji każdego z obrazów, tj. dla każdego próbki o współrzędnych <math>(x,y)</math>
-<math>IQ(x,y)=RF(x,y)\cdot e^{-i\cdot2\pi f_{0}t} </math><br><br>
+<math>IQ(x,y)=RF(x,y)\cdot e^{-i\cdot2\pi f_{0}t} </math>
 gdzie <math>IQ</math> - sygnał po demodulacji; <math>RF</math> - sygnał przed demodulacją; <math>t</math> - czas odpowiadający momentowi akwizycji próbki z danej głębokości; możemy przyjąć uproszczone założenie, że:
 <math>t=2y/c </math>; jak widać, pierwszy wymiar (szerokość) jest w naszej procedurze nieistotny.
@@ Linia 94: / Linia 95: @@
 ===Zadanie===
-Proszę przygotować funkcję estymującą częstotliwość średnią w oparciu o powyższy estymator dla jednowymiarowego zespolonego sygnału.
+Proszę przygotować funkcję estymującą częstotliwość średnią w oparciu o powyższy estymator dla jednowymiarowego sygnału zespolonego.
 ==Prezentacja Kolor==
 Po estymacji średnich prędkości w obszarach odpowiadających punktom na siatce użytej do rekonstrukcji obrazu, możemy nałożyć taką mapę na obraz B-mode w celu uzyskania obrazu "Kolor". Mając tablicę z danymi B-mode oraz tablicę prędkości, możemy otrzymać połączony obraz za pomocą poniższego skryptu:
 <source lang=python>
-#BMode - tablica zrekonstruowanych danych po przefiltrowaniu, obwiedni itp.
+# BMode - tablica zrekonstruowanych danych po przefiltrowaniu, obwiedni itp.
-#flow - tablica przepływów
+# flow - tablica przepływów
 Frame = Image.fromarray(np.uint8(cm.bone(BMode)*255))
 flowMask = np.copy(flow)
@@ Linia 115: / Linia 116: @@
 Frame.paste(flow, (0,0), flow)
 </source>
 ===Zadanie===
 Zaimplementować funkcję otrzymującą na wejściu trójwymiarową tablicę zawierającą obrazy RF z <math>N</math> kolejnych chwil pomiarowych i zwracającą obrazy z nałożoną na nią mapą prędkości obliczonych przy użyciu estymatora autokorelacyjnego.
@@ Linia 122: / Linia 124: @@
 ===Rozpoznawanie ruchu===
-Jedną z pierwszy obserwacji jakiej można dokonać, to zauważenie, że nasza mapa jest niezerowa w punktach w których nie spodziewamy się żadnego ruchu. Pewnym rozwiązaniem tego problemu może być zastosowanie prostego kryterium rozpoznawania ruchu. Możemy przyjąć, że sygnał pochodzący od struktur pozostających w spoczynku powinien być stały w czasie. W praktyce sygnały takie różnić będą się głównie o składową szumu elektronicznego. Jednocześnie, sygnał pochodzący od ruchomych struktur będzie charakteryzować się względnie dużą zmiennością w czasie. <br>
+Jedną z pierwszy obserwacji jakiej można dokonać, to zauważenie, że nasza mapa jest niezerowa w punktach w których nie spodziewamy się żadnego ruchu. Pewnym rozwiązaniem tego problemu może być zastosowanie prostego kryterium rozpoznawania ruchu. Możemy przyjąć, że sygnał pochodzący od struktur pozostających w spoczynku powinien być stały w czasie. W praktyce sygnały takie różnić będą się głównie o składową szumu elektronicznego. Jednocześnie, sygnał pochodzący od ruchomych struktur będzie charakteryzować się względnie dużą zmiennością w czasie.
 ====Zadanie====
 Zaimplementować procedurę zerującą estymatę prędkości w punkcie, jeśli średnia różnica między wartościami sygnału w tym punkcie w kolejnych chwilach czasu jest względnie mała (próg proszę dobrać eksperymentalnie - zaczynając np. od progu 10% średniej różnicy w obrazie).
 ===Filtr medianowy===
-Dobrym rozwiązaniem w tego typu przypadkach jest zastosowanie filtru medianowego. Filtr medianowy przekształca wartość w punkcie na medianę wartości sygnału w ustalonej liczbie sąsiednich próbek (w obu wymiarach):<br><br>
+Dobrym rozwiązaniem w tego typu przypadkach jest zastosowanie filtru medianowego. Filtr medianowy przekształca wartość w punkcie na medianę wartości sygnału w ustalonej liczbie sąsiednich próbek (w obu wymiarach):
-<math>s_{med}(x,y) = \sum^{K-1}_{i=0}\sum^{K-1}_{j=0}s(x+i,y+j) </math><br><br>
 gdzie <math>K</math> - rozmiar okna filtru.
 Filtr taki usuwa skrajne wartości w dużo większym stopniu niż np. filtr oparty o średnią arytmetyczną. W bibliotece scipy istnieje gotowa funkcja filtrująca tablicę filtrem medianowym:
 <source lang=python>
 scipy.signal.medfilt2d(array, K)
 </source>
 ====Zadanie====
-Porównać mapy prędkości przed i po filtracji medianowej (dla kilku różnych rozmiarów okna filtracji)
+Porównać mapy prędkości przed i po filtracji medianowej (dla kilku różnych rozmiarów okna filtracji).

@@ Linia 72: / Linia 72: @@
 i stąd
-<math> \dot{R}(0)=jA(0)\dot{\phi}(0) </math> oraz <math> R(0)=A(0) </math>.
+<math> \dot{R}(0)=jA(0)\dot{\phi}(0) </math>
 Korzystając ze wcześniejszych równań dostajemy
@@ Linia 82: / Linia 86: @@
 ====Estymator Millera-Rochwargera====
 Wartość <math>R(T_{PRF})</math> estymować możemy np. w oparciu o estymator Millera-Rochwargera <ref>Miller, Kenneth, and M. Rochwarger. "A covariance approach to spectral moment estimation." IEEE Transactions on Information Theory 18.5 (1972): 588-596.</ref>. Jest to [[wikipedia:Maximum_likelihood_estimation|estymator maksymalnej wiarygodności]] (maximum likehood estimator). Estymator ten jest skonstruowany dla par obserwacji zespolonego procesu <math>s</math> próbkowanych w równych odstępach czasu (w naszym przypadku <math>T_{PRF}</math>). W naszym przypadku jako pary obserwacji traktuje się pary próbek z kolejnych par następujących po sobie pomiarów (oddzielonych wartością <math>T_{PRF}</math>). Otrzymujemy wtedy oszacowanie wartości autokorelacji w oknie czasowym długości <math>N-1</math>
-<br><br><math>\bar{\omega}=\frac{1}{T_{PRF}}\mbox{arctan}\frac{Im(\sum^{N-2}_{i=0}s(i+1)\cdot \overline{s(i)})}{Re(\sum^{N-2}_{i=0}s(i+1)\cdot \overline{s(i)})}  </math><br><br>
 Pamiętając zależność między średnią częstotliwością <math>f_s</math> a średnią częstością kołową <math>fs=(1/2\pi)\bar{\omega}</math> dostajemy
-<br><br>
-<math>f_s=\frac{1}{2\pi T_{PRF}}\mbox{arctan}\frac{Im(\sum^{N-2}_{i=0}s(i+1)\cdot \overline{s(i)})}{Re(\sum^{N-2}_{i=0}s(i+1)\cdot \overline{s(i)})}  </math><br><br>
+<math>f_s=\frac{1}{2\pi T_{PRF}}\mbox{arctan}\frac{Im(\sum^{N-2}_{i=0}s(i+1)\cdot \overline{s(i)})}{Re(\sum^{N-2}_{i=0}s(i+1)\cdot \overline{s(i)})}  </math>
 ===Zadanie===

@@ Linia 36: / Linia 36: @@
 Dotychczas traktowaliśmy rekonstruowane obrazy jako sygnały dwuwymiarowe. Teraz, do zmierzenia średniej prędkości obiektów poruszających się w danym punkcie, konieczna będzie analiza wielu następujących po sobie obrazów, które zostały zebrane ze stałym interwałem czasowym (tzw. PRF - Pulse Repetition Frequency). W poniższych rozważaniach pomijamy wymiary przestrzenne - całe rozumowania przeprowadzone są dla sygnału <math>s(i)</math> na ustalonej głębokości i szerokości. Należy mieć na uwadze, że w praktyce ciężko mówić o estymacji prędkości w punkcie, ponieważ w naszej procedurze estymacji uwzględniamy efektywnie informacje z pewnego obszaru pomiarowego (energia danego piksela zawiera informację o obiektach z pewnego obszaru).
 Średnią prędkość będziemy mogli obliczyć korzystając ze średniego przesunięcia dopplerowskiego <math>\Delta f</math> w oparciu o szkolny wzór na [[Fizyka_III/Fale_dźwiękowe#Efekt_Dopplera|częstotliwość Dopplera]]:
 <math>\Delta f=\frac{2f_{0}v\mbox{cos}\theta}{c}</math>
 gdzie <math>v</math> średnia prędkość w punkcie pomiarowym; <math>\theta</math> kąt między kierunkiem przepływu a wiązką nadawczą.
 Przesunięcie dopplerowskie jest różnicą między średnią częstotliwością <math>f_0</math> sygnału nadanego a średnią częstotliwością <math>f_s</math> sygnału odebranego (pochodzącego od rozpraszaczy przemieszczających się w kierunku do lub od głowicy). Przypomnijmy, że zdemodulowany sygnał jest sygnałem o widmie pasmowym położonym w pobliżu częstości zerowej. Pozwala nam to przyjąć, że średnie przesunięcie częstotliwości odpowiada średniej częstotliwości tego sygnału <math>\Delta f = f_s-f_0 = f_s</math>.
@@ Linia 45: / Linia 47: @@
 Jednocześnie, przy założeniu słabej stacjonarności sygnału dopplerowskiego, na mocy [[Twierdzenie_Wienera-Chinczyna|twierdzenia Wienera-Chinczyna]], funkcję autokorelacji <math>R(\tau)</math> tego sygnału możemy wyrazić jako:
 <math> R(\tau)=\int^{\infty}_{-\infty} P(\omega)e^{j\omega\tau}\mbox{d}\omega</math>
 Różniczkując powyższe dostajemy:
 <math> \dot{R}(\tau)=j\int^{\infty}_{-\infty} \omega P(\omega)e^{j\omega\tau}\mbox{d}\omega</math>
 oraz
 <math> \ddot{R}(\tau)=-\int^{\infty}_{-\infty} \omega^{2} P(\omega)e^{j\omega\tau}\mbox{d}\omega</math>.
 Stąd, częstotliwość średnia możemy być również przedstawiona jako
 <math> \bar{\omega}=-j\frac{\dot{R}(0)}{R(0)}</math>
 W celu uproszczenia obliczeń przyjmuje się często następujące uproszczenie
 <math> R(\tau)= |R(\tau)|e^{j\phi(\tau)}=A(\tau)e^{j\phi(\tau)} </math>,
 gdzie <math>A(\tau)</math> jest funkcją parzystą i <math>e^{j\phi(\tau)}</math> jest funkcją nieparzystą. Wtedy
 <math>\dot{R}(\tau)=(\dot{A}(\tau)+jA(\tau)\dot{\phi}(\tau)e^{j\phi(\tau)} </math>
 i stąd
 <math> \dot{R}(0)=jA(0)\dot{\phi}(0) </math>
 oraz
 <math>R(0)=A(0) </math>
 Korzystając ze wcześniejszych równań dostajemy
 <math>\bar{\omega}=\dot{\phi}(0)\approx \frac{\phi(T_{PRF})-\phi(0)}{T_{PRF}}=\frac{\phi(T_{PRF})}{T_{PRF}}=\frac{1}{T_{PRF}}\frac{\mbox{Im}(R(T_{PRF}))}{\mbox{Re}(R(T_{PRF}))} </math>,
 gdzie <math>T_{PRF}</math> - czas między kolejnymi strzałami (odwrotność PRF).