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.