Skip to content

W czasopiśmie Otorynolaryngologia opublikowana została praca:

System automatycznej detekcji słuchowych potencjałów wywołanych pnia mózgu. I. Opis i testowanie systemu
Bartosz Trzaskowski, W. Wiktor Jędrzejczak, Edyta Piłka, Krzysztof Kochanek, Henryk Skarżyński
Otorynolaryngologia, 2013; 12(3): 137-147
http://www.mediton.pl/PL/czasopisma/otorynolaryngologia/archiwum.html

Streszczenie
Wprowadzenie. Słuchowe potencjały wywołane pnia mózgu (Auditory Brainstem Response, ABR) stanowią obecnie jedną z najczęściej stosowanych metod obiektywnego badania słuchu. Jednak subiektywna ocena wyniku badania wpływa na ryzyko wystąpienia różnic w ocenach prowadzonych przez różne osoby. Opracowanie systemu umożliwiającego wykrywanie odpowiedzi w sposób automatyczny mogłoby mieć pozytywny wpływ na standaryzację i obiektywizację procesu oceny badań ABR.

Cel pracy. Opracowanie systemu automatycznej detekcji słuchowych potencjałów wywołanych pnia mózgu oraz ocena jego działania.

Materiał i metody. System został oparty na metodach z zakresu analizy sygnałów i statystyki. Jego działanie zostało przetestowane na dużej liczbie danych symulowanych oraz na rzeczywistych sygnałach zarejestrowanych w grupie osób z prawidłowym słuchem oraz w grupie pacjentów z ubytkami słuchu.

Wyniki. Różnice w oznaczeniach systemu względem oznaczeń ekspertów nie były znacząco wyższe od typowo spotykanych w praktyce klinicznej różnic w oznaczeniach pomiędzy ekspertami.

Wnioski. Otrzymane wyniki sugerują, że prezentowany algorytm może być użyteczny jako system wspomagający pracę osoby oceniającej badanie ABR.

W poprzednim wpisie przedstawiłem funkcję odpowiedzi aparatu Nikon D70 przy zapisie zdjęć w formacie JPG Fine. Tutaj chciałbym pokazać jak wygląda funkcja odpowiedzi w przypadku zapisu obrazu w formacie NEF, będącym natywnym formatem Nikona do zapisu danych surowych (RAW). Dane w tym formacie są zasadniczo nieprzetworzonym zapisem tego co zarejestrowała matryca aparatu. Pisząc nieprzetworzonym mam na myśli w tym przypadku brak przetworzenia w postaci konwersji gamma, czy balansu bieli. Jednak przetworzenie obrazu rejestrowanego na matrycy musi mieć miejsce już na etapie jego digitalizacji czyli przekształcenia informacji analogowej z matrycy (w postaci ładunku) na informację w postaci cyfrowej. Taka kwantyzacja ciągłej wartości wiąże się ze zmianą informacji. Ale jak donoszą niektórzy dociekliwi użytkownicy, Nikon miesza już nawet na tym etapie przetwarzając i zmieniając informację o obrazie.

W szczególności udowodniono, że zapis zdjęć w formacie NEF z bezstratną kompresją stosowany w aparacie Nikon D70, jest wbrew podawanym przez producenta informacjom stratny. W głównej mierze, strata informacji związana jest z obniżeniem rozdzielczości w jasnej części zakresu. Nikon D70 jest wyposażony w matrycę Sony ICX413AQ i 12-bitowy przetwornik analogowo-cyfrowy. Rozdzielczość 12-bitowa pozwala na zapis 2^12=4096 poziomów jasności, ale podczas konwersji do formatu RAW liczba poziomów jest limitowana do 683 i dopiero potem poddawana bezstratnej kompresji słownikowej podobnej do stosowanej w plikach ZIP. O ile kompresja rzeczywiście jest bezstratna to jednak informacja tracona jest na etapie kwantyzacji do 683 dyskretnych wartości jasności. Krzywa kwantyzacji jest zapisywana w plikach NEF. Kodowanie to zachowuje pełen zakres jasności, jednak zapis 12-bitowej informacji (4096 poziomów) w 683 dyskretnych wartościach powoduje zmniejszenie rozdzielczości w jasności. Kształt krzywej kwantyzacji, na początku liniowy, a następnie wzrastający z kwadratem, powoduje że rozdzielczość maleje wraz ze wzrostem jasności. Zastosowanie tego rodzaju konwersji miało prawdopodobnie na celu znaczne przyspieszenie (prawie o rząd wielkości) zapisu plików NEF przez aparat. Starsze modele (jak np D1H czy D100) potrafiły przetwarzać obraz przez 20-30 sekund wykonując kompresję przed zapisem. Obecne część nowych modeli aparatów posiada możliwość wyboru trybu zapisu plików NEF pomiędzy 12- i 14-bitowymi, oraz skompresowanymi i nieskompresowanymi.

Dodatkowo pewna utrata informacji w D70 związana jest z zastosowaniem optycznego filtra dolnoprzepustowego, mającego na celu usunięcie z obrazu składowych wysokoczęstotliwościowych. Matryca aparatu podczas zapisu obrazu niejako próbkuje scenę z pewną określoną częstotliwością zwaną częstotliwością próbkowania. Jeżeli w fotografowanej scenie występują jednak częstotliwości większe od połowy częstotliwości próbkowania to w rejestrowanym obrazie wystąpi zjawisko aliasingu, fotografom bądź grafikom znane pod nazwą efektu mory. Żeby przeciwdziałać temu zjawisku, producenci przez wiele lat umieszczali przed matrycą filtr rozmywający (usuwający wysokie częstotliwości) dostosowany do danej matrycy. To rozwiązanie gwarantowało, że częstotliwości wyższe niż dopuszczalne dla danej matrycy będą odfiltrowywane, a matryca zarejestruje tylko dozwolone częstotliwości. Łatwo też jest się domyślić jaka jest zależność częstotliwości próbkowania od rozmiaru matrycy (w megapikselach). Im większy rozmiar w MPx, tym większa jest częstotliwość próbkowania, czyli tym wyższe częstotliwości mogą być przez daną matrycę zarejestrowane bez ryzyka wystąpienia efektu aliasingu. Ostatnio, w związku postępem technologicznym i ciągłym wzrostem rozmiarów rejestrowanych obrazów, daje się zauważyć tendencję do usuwania tego filtru co bardzo pozytywnie wpływa na poprawę ostrości i zwiększenie szczegółowości zdjęć. Również Nikon w swoim najnowszym modelu D5300 zrezygnował z optycznego filtru dolnoprzepustowego uznając, że problem efektu mory nie jest przy rozmiarze zdjęcia 24 MPx już tak bardzo istotny.

Niemniej jednak, pomimo wspomnianych wyżej przekształceń, można by się spodziewać, że ponieważ w aparacie nie została wykonana kompresja gamma (jak w przypadku plików JPG), a matryca w Nikonie D70 to matryca CCD (dla której ładunek jest proporcjonalny do ekspozycji) to zależność odpowiedzi może być tutaj zbliżona do liniowej. I to przede wszystkim chciałem sprawdzić wykonując przedstawione poniżej obliczenia.

Przeanalizowałem dwie serie zdjęć wykonanych z różnymi ustawieniami ekspozycji. Zdjęcia były wykonane dokładnie tak samo jak opisałem w poprzednim artykule, z tą różnicą, że tym razem zapisane zostały w formacie NEF. Są to dokładnie te same sceny, sfotografowane w tych samych warunkach, ponieważ serie zdjęć NEF były wykonywane bezpośrednio przed seriami JPG.

Tak wyglądają serie zdjęć w formacie NEF:

Barckets. Series 1. NEF.

Barckets. Series 1. NEF.

Barckets. Series 2. NEF.

Barckets. Series 2. NEF.

Łatwo zauważyć, że w przypadku obydwu serii, zdjęcia w formacie NEF wyglądają na dużo ciemniejsze i o większym kontraście niż zdjęcia w formacie JPG z poprzedniego artykułu. Dzieje się tak dlatego, że zdjęcia w formacie NEF są w liniowej przestrzeni kolorów, a zdjęcia w formacie JPG w logarytmicznej. Nasz narząd wzroku podobnie jak inne nasze zmysły działa w sposób naturalny logarytmicznie. Zdjęcia w logarytmicznej przestrzeni kolorów (np wykonane przy pomocy tradycyjnej fotografii wykorzystującej kliszę fotograficzną) wyglądają dla nas naturalnie i fotorealistycznie. Dlatego aparat przed zapisem zdjęć w docelowym formacie JPG skompresował je do logarytmicznej przestrzeni kolorów wykonując przekształcenie potęgowania – tzw konwersję gamma. To spowodowało bardziej naturalnie dla nas wyglądający rozkład jasności pikseli w zdjęciach.

Faktyczna rozdzielczość bitowa w plikach NEF zapisywanych w aparacie Nikon D70 wynosi: log2(683)=9.42 bity. Żeby wykonać obliczenia przekształciłem pliki NEF w 8-bitowe pliki TIFF. W tym celu wykorzystałem napisany przez Dave Coffina program dcraw. Konwersja nie mogła wykonywać żadnych automatycznych zmian jasności (domyślnie dcraw rozciąga histogram w ten sposób, żeby 1% pikseli był wyświetlany jako biały), wykonywać konwersji gamma, repróbkować pikseli, ani ingerować w przestrzeń kolorów. Polecenie którego użyłem to: „dcraw -T -W -g 1 1 -v -j -o 0″. Zdjęcia z Nikona D70 w formacie NEF mają rozmiar 3039×2014 pikseli. Maska pikseli wybranych do obliczeń była utworzona identycznie jak w przypadku obrazów JPG, czyli tablica 30×20 pikseli równomiernie rozmieszczonych w obrazie z zachowaniem 5% marginesu od krawędzi obrazu.

Dygresja:
Żeby poprawnie oszacować przy pomocy tej metody krzywą odpowiedzi dla np 14-bitowych plików NEF, dla serii 11 ekspozycji, zgodnie z zależnością: N*(P-1)>(Zmax-Zmin), trzeba by było wykorzystać w obliczeniach: N>(2^14-1)/(11-1) czyli przynajmniej 1639 pikseli. Tablica układu równań liniowych zajmowałaby wówczas w pamięci: (1639*11+2^14+1)*(2^14+1639)*16/1000/1000/1000 = 9.9 GB RAM. Oczywiście funkcja zaproponowana przez Debeveca i Malika musiałaby zostać zmodyfikowana w celu uwzględnienia 14-bitowej rozdzielczości: Zmax=2^14-1 i przesunięcia środka rozkładu do jasności Z=(2^14)/2-1.
Koniec dygresji.

Krzywe odpowiedzi dla plików NEF z aparatu Nikon D70:

Nikon D70 response curve. NEF

Nikon D70 response curve. NEF

Nikon D70 response curve. NEF

Nikon D70 response curve. NEF

Widać, że otrzymana krzywa odpowiedzi dla zapisu NEF ma bardziej stromy przebieg niż krzywa dla kodowania JPG z poprzedniego artykułu. Jednak trudno jest coś więcej powiedzieć, ponieważ na tych ilustracjach wynik jest przedstawiony na płaszczyźnie pół-logarytmicznej. Żeby dokładniej porównać obie funkcje mamy dwie możliwości: pokazać je na wykresie logarytmiczno-logarytmicznym lub na wykresie liniowo-liniowym.

Wykres o obu osiach liniowych dla serii 1:

Nikon D70 response curve. NEF. lin-lin

Nikon D70 response curve. NEF. lin-lin

Wykres dla obu osi logarytmicznych dla serii 1:

Nikon D70 response curve. NEF. log-log

Nikon D70 response curve. NEF. log-log

Na obu tych wykresach widać, że w przypadku zapisu w formacie NEF zależność aż do momentu wysycenia jest rzeczywiście liniowa.

Na wykresie log-log wątpliwości co do jakości wyniku oszacowania krzywej odpowiedzi, mogą budzić punkty zlokalizowane w lewej dolnej części wykresu. Duży rozrzut punktów wzdłuż osi X (logarytmiczna wartość ekspozycji) oraz duże odstępy pomiędzy kolejnymi wartościami na osi Y (logarytmiczna wartość jasności pikseli). Kwestia dużych odstępów wzdłuż osi Y stanie się jasna gdy tylko zdamy sobie sprawę co dokładnie powoduje transformacja skali liniowej do logarytmicznej. Skala logarytmiczna sprawia, że niskie wartości są na osi rozciągane, a wyższe coraz bardziej kompresowane. Czyli duże duże odstępy w tym rejonie to naturalna cecha przekształcenia wartości w skali liniowej do logarytmicznej. Możemy policzyć, że do połowy tego wykresu czyli do wartości 3 na osi Y, znajduje się exp(3)=20 pikseli spośród wszystkich 256. Czyli w tej skali, połowę wykresu zajmuje niecałe 8% wszystkich możliwych wartości pikseli. W skali liniowej te wartości zajmowałyby znikomą powierzchnię w lewym dolnym rogu płaszczyzny wykresu. A kwestię rozrzutu wzdłuż osi X można wyjaśnić uwzględniając fakt, że matryce CCD posiadają naturalną tendencję do rejestrowania szumu przy małych wartości ekspozycji, a połowę wykresu zajmują właśnie piksele dla najniższych 20 wartości ekspozycji. Dlatego też w procesie estymacji krzywej odpowiedzi stosowana był funkcja zmniejszająca wagę tych wartości.

Ale jak wygląda na tych wykresach zależność krzywej odpowiedzi dla zapisów JPG?
Wykres lin-lin:

Nikon D70 response curve. JPG. lin-lin

Nikon D70 response curve. JPG. lin-lin

Wykres log-log:

Nikon D70 response curve. JPG. log-log

Nikon D70 response curve. JPG. log-log

Na obu wykresach dla formatu JPEG widać dla krzywych każdego koloru charakterystyczną nieliniową zależność przypominającą transformację gamma.

Z powyższych dwóch artykułów i przedstawionych wyników można wyciągnąć wnioski dotyczące aparatu Nikon D70:

  • zapis obrazów w formacie JPEG, wiąże się z zastosowaniem nieliniowego przekształcenia o profilu pokazanym na powyższych ilustracjach
  • zapis w formacie NEF i wykonanie płaskiej konwersji do 8-bitowych obrazów TIFF w programie dcraw, daje obrazy z liniową (w pewnym zakresie) funkcją odpowiedzi.

W 1997 roku Debevec i Malik opublikowali interesującą pracę. Przedstawili w niej między innymi metodę pozwalającą na oszacowanie funkcji odpowiedzi systemu powstawania obrazu w oparciu o serię zdjęć wykonanych przy różnych wartościach ekspozycji. Metoda ta polega na znalezieniu funkcji minimalizującej w sensie metody najmniejszych kwadratów błąd rozwiązania układu równań liniowych wiążących wartość jasności pikseli obrazu z ekspozycją.

W celu określenia funkcji odpowiedzi cyfrowej lustrzanki Nikon D70 dla zdjęć zapisywanych przez aparat w formacie JPG wykonałem trzy serie ujęć. Każda seria składała się z jedenastu zdjęć wykonywanych dla kolejnych ustawień ekspozycji z krokiem 1EV w przedziale [-5EV, 5EV] względem właściwej wartości nastawu.

Do zdjęć zostały wybrane sceny nie zawierające elementów ruchomych i neutralne kolorystycznie. Cały kadr zawierał obiekty o szarych kolorach bez elementów o dużej saturacji. Fotografowane sceny zawierały rejony kontrastujące w jasności aby zapewnić szeroki rozkład punktów w przestrzeni [ekspozycja, jasność piksela] dla pojedynczych wartości ekspozycji. Miało to na celu poprawę jakości dopasowywania do siebie krzywych dla pojedynczych pikseli w trakcie obliczeń metodą najmniejszych kwadratów.

Przykładowe dwie serie zdjęć użyte do obliczeń przedstawione są na rysunkach poniżej.

Seria 1:

Brackets. Series 1. JPG

Brackets. Series 1. JPG

Seria 2:

Brackets. Series 2. JPG

Brackets. Series 2. JPG

Zdjęcia były wykonywane przy użyciu statywu, kolejno jedno po drugim, w odstępach czasu najkrótszych na jakie pozwalał aparat kontrolowany za pomocą kabla USB OTG przez program zdalnie wyzwalający migawkę zainstalowany na smartfonie. Wykonanie zdjęć w możliwie najkrótszym czasie było ważne żeby zapobiec zmianie warunków oświetleniowych podczas fotografowania. Wszystkie serie zdjęć zostały zrobione w pochmurny dzień przy słońcu całkowicie zasłoniętym chmurami. Zdjęcia wykonywano zachowując stałą liczbę przesłony, uzyskując zmienną wartość ekspozycji operując czasem otwarcia migawki. Pozwoliło to uniknąć problemów związanych ze zmianą głębi ostrości i winietowania. W celu zmniejszenia wkładu szumu z matrycy w rejonach o niskiej jasności, zdjęcia wykonywano przy najniższej możliwej wartości czułości 200 ISO. Wielkość zdjęć ustawiono na maksymalną.

Jako czasy otwarcia migawki uwzględnione w obliczeniach użyte zostały wartości wyświetlane przez aparat czyli np 1/500, 1/250 czy 1/125 s, pomimo tego, że Debevec i Malik sugerowali, że czasy będące wynikiem potęgowania liczby 2 (odpowiednio: 1/512, 1/256 czy 1/128 s) są bardziej zgodne z rzeczywistym czasem ekspozycji.

Z punktu widzenia stosunku jakości do wydajności obliczeń istotny jest wybór liczby pikseli obrazu dla których obliczenia będą wykonywane. Nikon D70 posiada 6 MPx matrycą, a największe zdjęcia zapisywane przez aparat w jakości JPG Fine są w rozmiarze 3008×2000 px. Ponieważ układ równań liniowych minimalizujących błąd znalezienia funkcji odpowiedzi jest rzędu N*P + Zmax-Zmin to wykorzystanie wszystkich 6.016.000 pikseli, przy serii 11 ekspozycji, dla obrazów 8-bit, wymagałoby utworzenia w pamięci tablicy o rozmiarze ok (3008*2000*11+2^8+1)*(2^8+3008*2000)*8/1000/1000/1000=3,185,066 GB pamięci RAM. Być może w przyszłości będzie można wykonywać takie obliczenia na telefonie komórkowym jednak obecnie nie są one (wg wiedzy autora) możliwe do przeprowadzenia. W rzeczywistości, żeby określić ten układ równań wystarczy uwzględnić N pikseli spełniających nierówność: N*(P-1)>(Zmax-Zmin); gdzie N – oznacza liczbę pikseli; P – liczbę ekspozycji; Zmax-Zmin to maksymalna różnica w jasności pikseli. Przy 11 ekspozycjach dla obrazu 8-bitowego, wystarczającą liczbą pikseli do obliczeń będzie: N>(2^8-1)/(11-1) czyli N>=26.

Znając wymaganą minimalna liczbę pikseli do obliczenia krzywej odpowiedzi, pozostaje jeszcze określenie sposobu próbkowania obrazu. W pracy z 1997 roku Debevec i Malik ręcznie określali które piksele miały być uwzględnione w obliczeniach. Ja zdecydowałem się na wybranie do obliczeń tablicy 600 pikseli (30×20) rozmieszczonych w jednorodnych odstępach od siebie z zachowaniem 5% marginesu od krawędzi obrazu. Ze względu na to że liczba pikseli była o rząd wielkości większa niż wymagana, układ charakteryzował się dużą redundancją, a w pamięci zajmował:
(30*20*11+2^8+1)*(2^8+30*20)*8/1000/1000=46.96MB RAM. Obliczenia na typowych komputerach były kwestią sekund.

Poniżej pokazane są wykresy otrzymanych funkcji odpowiedzi Nikona D70 dla zapisu JPG Fine dla trzech kolorów. W celu zmniejszenia wpływu skrajnych wartości jasności na oszacowanie funkcji, zastosowane zostało okno trójkątne jako funkcja określająca wagę jasności pikseli. Na ilustracjach przedstawione są także punkty w przestrzeni [ekspozycja, wartość piksela] które posłużyły do estymacji funkcji odpowiedzi.

Krzywa odpowiedzi aparatu Nikon D70 przy zapisie w formacie JPEG dla serii 1:

Nikon D70 response curve. JPG

Nikon D70 response curve. JPG

Krzywa odpowiedzi aparatu Nikon D70 przy zapisie w formacie JPEG dla serii 2:

Nikon D70 response curve. JPG

Nikon D70 response curve. JPG

Otrzymane krzywe są praktycznie identyczne. Przedstawiają poszukiwaną funkcję odpowiedzi aparatu.

Strona została przeniesiona na nowy serwer. Wszystkie błędy proszę zgłaszać tutaj.

Na nowej stronie Tekstury dostępne są do pobrania darmowe tekstury. Tekstury wysokiej jakości dostępne są odpłatnie.

W czasopiśmie Computers in Biology and Medicine ukazał się artykuł opisujący system automatycznej detekcji i usuwania fali sonomotorycznej z zapisów słuchowych potencjałów wywołanych pnia mózgu (ABR): Automatic removal of sonomotor waves from auditory brainstem responses

Streszczenie:

We have developed a computerized technique for automatic detection and removal of sonomotor waves (SMWs) from auditory brainstem responses (ABRs). Our approach is based on adaptive decomposition using a redundant set of Gaussian and 1-cycle-limited Gabor functions. In order to find optimal parameters and evaluate the efficiency of the methods, simulated data were first used before applying it to clinical data. Results were good and confirmed by an expert with years of clinical experience in ABR evaluation.

Z okazji zbliżających się Świąt Bożego Narodzenia i Nowego Roku 2013 serdeczne życzenia wszystkiego dobrego dla wszystkich odwiedzających tę stronę.

Bartosz Trzaskowski

W tym wpisie chciałbym porównać podstawowe cechy dwóch platform do obliczeń naukowych: środowiska Matlab firmy Mathworks i języka Python z dedykowanymi do tego typu zadań bibliotekami. Prosty program obliczający widmo sygnału powinien być odpowiedni do porównania podstawowych cech obu języków, podobieństw i różnic składni, łatwości implementacji, przejrzystości i zrozumiałości kodu.

Poniższy skrypt w języku Matlaba liczy widmową gęstość mocy (PSD) modelowanego sygnału zawierającego dwie składowe harmoniczne o częstotliwościach 30 i 80 Hz. Widmo liczone jest przy wykorzystaniu szybkiej transformaty Fouriera, okienkowanie oknem Hanninga.

tic
 
fs=400; T=10;
 
nfft=512;
overlap=0.5;
 
A1=1; f1=30; fi1=0;
A2=2; f2=80; fi2=2;
An=1.5;
 
t=0:1/fs:T;
t=t(:); % make 't' a column vector
 
sig=A1*cos(2*pi*f1*t+fi1) + A2*cos(2*pi*f2*t+fi2);
noise=An*rand(length(t),1);
 
x=sig+noise;
 
wind=hanning(nfft);
 
noverlap=fix(nfft*overlap);
nadvance=nfft-noverlap;
nrecs=fix((length(t)-noverlap)/nadvance);
 
Pxx=zeros(nfft,1);
ind=1:nfft;
for k=1:nrecs
    xs=x(ind);
    xs=(xs(:)-mean(xs)).*wind;
    Xf=fft(xs,nfft)/nfft;
    Pxx=Pxx+abs(Xf).^2;
    ind=ind+nadvance;
end
 
Pxx=Pxx(1:nfft/2);
f=(0:nfft/2-1)/(nfft/2-1)*(fs/2);
 
figure
plot(f,10*log10(Pxx)) % in dB
xlabel('f [Hz]')
ylabel('PSD [dB]')
saveas(gcf,'fig_PSD_matlab.png','png')
 
% time
czas=toc;
czas_=['Time: ' int2str(fix(czas/3600)) ' hours, ' ...
    int2str(fix(czas/60)-fix(czas/3600)*60) ' min and ' ...
    num2str(mod(czas,60),'%0.2f') ' s.'];
disp(' ');
disp(czas_)

Ten sam program napisany w Pythonie przy wykorzystaniu bibliotek NumPy i Matplotlib wygląda tak:

import numpy, pylab, time
from math import pi
 
tic=time.clock()
 
fs,T=400,10 
 
nfft=512
overlap=0.5
 
A1,f1,fi1 = 1,30,0 
A2,f2,fi2 = 2,80,2 
An=1.5
 
t=numpy.arange(0,T,1./fs)
t=t.transpose()
 
sig=A1*numpy.cos(2*pi*f1*t+fi1) + A2*numpy.cos(2*pi*f2*t+fi2)
noise=An*numpy.random.rand(len(t))
x=sig+noise
 
wind=numpy.hanning(nfft)
 
noverlap=numpy.fix(nfft * overlap)
nadvance=nfft-noverlap
nrecs=numpy.fix((len(t)-noverlap)/nadvance)
 
Pxx=numpy.zeros((nfft))
ind=numpy.arange(0,nfft,1)
for k in range(nrecs):
    xs=x[ind]
    xs=numpy.multiply((xs-numpy.mean(xs)),wind)
    Xf=numpy.fft.fft(xs,nfft)/nfft
    Pxx=Pxx+numpy.square(numpy.abs(Xf))
    ind=ind+int(nadvance)
 
Pxx=Pxx[0:nfft/2]
f=numpy.arange(0.0,nfft/2,1)/(nfft/2-1)*(fs/2)
 
pylab.figure
pylab.plot(f,10*numpy.log10(Pxx),linewidth=0.5) # in dB
pylab.xlabel('f [Hz]')
pylab.ylabel('PSD [dB]')
pylab.savefig("fig_PSD_python.png")  
pylab.show()
 
# time
toc = time.clock()
czas=toc - tic
czas_='Time: '
czas_=czas_ + "%0.*f" %(0,numpy.fix(czas/3600)) + ' hours, '
czas_=czas_ + "%0.*f" %(0,numpy.fix(czas/60)-numpy.fix(czas/3600)*60) + ' min and '
czas_=czas_ + "%0.*f" %(2,numpy.mod(czas,60)) + ' s.'
print(czas_)

Na pierwszy rzut oka widać duże podobieństwo pomiędzy obydwoma programami i wydaje się, że dla użytkownika Matlaba rozpoczęcie pracy z Pythonem przy wykorzystaniu dedykowanych bibliotek powinno być względnie proste. Oba języki są językami wysokiego poziomu, dającymi dużą swobodę w wyborze metody rozwiązania problemu. Biblioteki do obliczeń numerycznych zawierają funkcje o bardzo podobnych nazwach, wywołaniu i działaniu do funkcji znanych użytkownikom Matlaba. Niemniej jednak, po kilku latach używania Matlaba, jest kilka elementów do których się przyzwyczaiłem i których jest mi brak przy pierwszym kontakcie z Pythonem. Na przykład: brak operatorów z kropką wyraźnie określających działania na pojedynczych elementach tablicy, czy przejrzystego i prostego zapisu macierzy w postaci a=[1,2,3;4,5,6].

Poniższy rysunek pokazuje wynik działania funkcji plot() w Matlabie (po lewej) i jej odpowiednika o takiej samej nazwie z biblioteki Matplotlib (po prawej).

Orientacyjny czas wykonania obydwu programów przy kilkukrotnym uruchomieniu jest bardzo zbliżony do siebie. Na przenośnym komputerze z procesorem Intel Core 2 Duo 2GHz wynosi:
Matlab: 0.25 s [0.24-0.26]
Python: 0.26 s [0.25-0.27]

Ciekawostką jest, że wykonanie skryptu Matlaba w środowisku Octave trwało znacznie krócej:
Octave: 0.08 s [0.07-0.09]

Żeby porównać wydajność obu środowisk, należałoby wykonać dużo większą ilość obliczeń. W następnym wpisie przedstawię wyniki takiego zestawienia wykorzystując programy wymagające znacznie większej mocy obliczeniowej.

Ostatnio coraz więcej osób pracujących naukowo, do wykonywania obliczeń numerycznych zaczyna wykorzystywać język Python . Ja również od pewnego czasu zastanawiam się nad przeniesieniem części swojej pracy z Matlaba do Pythona. Na korzyść Pythona przemawia kilka argumentów.

Przede wszystkim jest darmowy. W zestawieniu z kosztem Matlaba jest to istotny argument ponieważ znacząca jest już cena samego środowiska obliczeniowego, a dokupienie opcjonalnych toolboxów dodatkowo zwiększa koszt zakupu platformy Mathworks. Darmowe biblioteki Pythona jak np. NumPy, SciPy, Matplotlib czy IPython zapewniają dużą część funkcjonalności podstawowej wersji Matlaba i niektórych toolboxów.

Po drugie Python jest językiem programowania o interpreterach dostępnych na prawie wszystkie platformy. W wielu systemach, wraz z podstawowymi bibliotekami, jest domyślnie rozprowadzany razem z dystrybucją. Dzięki temu, po zainstalowaniu wymaganych bibliotek, możliwe jest wykonanie programów napisanych w Pythonie na prawie każdym komputerze. Natomiast środowisko Matlab, głównie ze względu na swój koszt, jest kupowane najczęściej przez instytucje naukowe, uczelnie lub duże firmy. I przeważnie z licencją na ograniczoną, niewielką liczbę stanowisk lub równoczesnych sesji. Dlatego w Matlabie, programy uruchomić będą mogły tylko osoby mające dostęp do tego środowiska, najczęściej w pracy lub na uczelni. W tym miejscu należałoby wspomnieć o darmowej alternatywie dla programów napisanych w języku Matlaba – platformie Octave. Język Octave jest prawie całkowicie zgodny z językiem Matlaba, a wszelkie różnice w składni i działaniu poszczególnych funkcji są uważane za błędy i systematycznie likwidowane. Jednak obie platformy nie są jeszcze w 100% zgodne ze sobą [Porting programs from Matlab to Octave] i o ile prostsze programy Matlaba bez problemu dają się uruchomić w Octave, to z wykonaniem bardziej skomplikowanych skryptów mogą wystąpić trudności. Umieszczone na tej stronie, proste skrypty w języku Matlaba powinny uruchamiać się w środowisku Octave w wersji 3.2.4 lub wyższej. Natomiast programy w języku Python, każdy powinien móc uruchomić na swoim komputerze (ewentualnie, po uprzednim zainstalowaniu Pythona).

Kolejnym argumentem na korzyść Pythona jest to, że jest on pełnym językiem programowania, a nie środowiskiem dedykowanym do przeprowadzania obliczeń. To powoduje, że w razie potrzeby łatwiej jest wykorzystać Pythona do wykonywania operacji nie związanych wprost z obliczeniami. Dodatkowo duża liczba dedykowanych bibliotek, optymalizowanych pod kątem realizacji określonych zadań bardzo upraszcza implementację. Python może być używany jako język skryptowy służący np. automatyzacji powtarzających się zadań, a jednocześnie po zaimportowaniu odpowiednich bibliotek można go w sposób prosty wykorzystać do przeprowadzania naukowych obliczeń. Wiele aplikacji udostępnia swoje API właśnie w Pythonie. Jestem w stanie bez większych problemów wyobrazić sobie sytuację, gdzie w celu wykorzystania części swoich wcześniejszych programów w jednej z takich aplikacji, przekopiowuję bez żadnych zmian funkcje i fragmenty programów napisanych w Pythonie. Skrypty stworzone w Matlabie są w takiej sytuacji bezużyteczne.

Dla równowagi podkreślić należy kilka zalet środowiska Mathworks. Matlab jest jedną z najszerzej wykorzystywanych platform obliczeniowych na świecie. Od wielu lat wykorzystywany jest przez środowiska naukowe w różnych dziedzinach badań i do dziś stanowi standard w środowiskach uniwersyteckich. Przy wykorzystaniu Matlaba, wykonywane są obliczenia do większości publikacji naukowych, a przykładowe fragmenty kodu i algorytmy zapisywane są w języku Matlaba. Chociażby z tego powodu, aby pozostać w zgodzie z tym standardem, warto jest przynajmniej znać zapis Matlabowy.

Inną zaletą platformy Mathworks względem Pythona jest to, że dostarcza jednorodne, kompletne proste w instalacji środowisko zawierające w jednym miejscu wszystkie wymagane do pracy elementy. Bezpośrednio po instalacji, podczas której automatycznie instalowane są również zakupione toolboxy, użytkownik uruchamia jedno zwarte środowisko w którym dostępne są: okno poleceń, zintegrowany edytor, podgląd zmiennych w przestrzeniach pamięci, debugger, profiler, historia poleceń, przeglądarka plików, pomoc on-line, a nawet dostęp do internetowych baz skryptów udostępnianych przez społeczność użytkowników Matlaba. Pracę ułatwiają między innymi kreatory importu i eksportu danych z popularnych formatów plików, obsługa schowka, czy zintegrowane narzędzia do edycji wykresów. Moim zdaniem, to jest bardzo duży atut Matlaba, że zaraz po krótkiej instalacji stanowi kompletne środowisko zawierające zintegrowane liczne narzędzia ułatwiające pracę.

Innym bardzo dużym atutem produktu Mathworks jest Simulink - środowisko do przeprowadzania symulacji i projektowania, zintegrowane z Matlabem. Do Simulinka również można dokupić dedykowane toolboxy. Simulink w połączeniu z Matlabem stanowi potężną platformę mogącą służyć analizie zachowań systemów dynamicznych, wizualizacji symulacji, modelowaniu danych i wielu innym zastosowaniom. Obecnie nie ma na rynku konkurencyjnego do Simulinka oprogramowania o podobnych możliwościach.

W internecie znajduje się wiele stron porównujących Matlaba i Pythona przy różnego rodzaju zastosowaniach. Wymienione powyżej argumenty za każdym z tych dwóch środowisk to moje subiektywne zdanie. Wypunktowałem różnice, które są istotne według mnie i w konkretnych zastosowaniach, z którymi mam do czynienia – obliczeniami w zakresie analizy sygnałów. Inne osoby mogą przywiązywać wagę do innych cech obu języków w zależności od swoich przyzwyczajeń lub oczekiwań. Żeby móc ocenić praktyczne aspekty pracy z tymi językami, chciałbym wykonać kilka obliczeń w obu środowiskach porównując łatwość implementacji, czas wykonania itp. Ponieważ z jednej strony nie mam dużego doświadczenia z programowaniem w Pythonie, a z drugiej strony chciałbym zrobić coś ciekawszego niż wypisanie na ekranie liczb w pętli, to w następnych artykułach przedstawię zastosowanie obu języków do kilku prostych przypadków analizy modelowanych sygnałów, a także dźwięku i obrazu.

Użytkownicy systemów Windows, mogą pobrać skompilowane wersje wielu bibliotek Pythona do obliczeń naukowych ze strony: Unofficial Windows Binaries for Python Extension Packages, utrzymywanej przez Christopha Gohlke. Użytkownikom systemów linuksowych najłatwiej będzie użyć bibliotek z oficjalnych repozytoriów lub skompilować najnowsze wersje ze źródeł.

Na zakończenie odnośniki do dwóch środowisk obliczeniowych o otwartym kodzie źródłowym: Scilab i SAGE.

W czasopiśmie Acta Bio-Optica et Informatica Medica Inżynieria Biomedyczna po dość długim oczekiwaniu ukazała się praca: Zastosowanie systemu automatycznej detekcji i usuwania fali sonomotorycznej z zapisu słuchowych potencjałów wywołanych pnia mózgu.

Publikacja opisuje wczesną wersję systemu. Praca opisująca bardziej rozbudowaną wersję algorytmu została złożona do druku w innym międzynarodowym czasopiśmie i w chwili obecnej oczekuje na recenzję.