W poprzednim wpisie omówiłem w jaki sposób generować liczby pseudolosowe w MATLABie. Skoro już wiemy jak to robić, czas zająć się analizą sygnałów stochastycznych w celu określenia ich właściwości. Dzisiaj na warsztat bierzemy MATLAB histogram – często wyznaczaną charakterystykę statystyczną, która daje nam dużo cennych informacji o sygnale!

Z tego posta wiemy, że w przypadku sygnałów losowych nie da się ściśle określić ich wartości w danej chwili czasu. Można jedynie oszacować częstość lub prawdopodobieństwo wystąpienia określonej wartości. Stąd narzędzia używane do analizy sygnałów losowych są specyficzne i nieco odmienne od tych, których używamy w przypadku sygnałów zdeterminowanych, na przykład sinusoidalnych.

MATLAB Histogram - podstawowe użycie

Podstawową charakterystyką określającą tzw. rozkład wartości sygnału jest histogram. Ilustruje on ile razy dana wartość wystąpiła w obserwowanym sygnale, a ściśle rzecz biorąc ile razy wystąpiły wartości z określonego przedziału. Na przykład: histogram sygnału stałego w czasie zawierałby tylko jeden prążek, odpowiadający tej stałej wartości; histogram sygnału binarnego posiadałby dwa prążki o wartości 0 i 1, i tak dalej. To na ile przedziałów dzielimy zakres zmienności sygnału zależy od użytkownika, a liczba tych przedziałów jest podstawowym parametrem histogramu nBins.

Do wyznaczenia histogramu w MATLABie przez wiele edycji służyła funkcja hist, która została zastąpiona od 2014 roku przez funkcję MATLABhistogram. Jej podstawowa składnia jest następująca:

historam(X, nBins);

gdzie X – to analizowany sygnał, a nBins – to liczba przedziałów. Od liczby przedziałów zależy czytelność i „szczegółowość” rysowanego histogramu.

Przeanalizujcie poniższy program. W zależności od wyboru rodzaju rozkładu (zmienna typ) generowany jest sygnał losowy o rozkładzie normalnym lub Chi-kwadrat i jest rysowany jego histogram. Proces jest powtarzany w pętli, a rysunek odświeżany jest  co sekundę. Za liczbę przedziałów odpowiada zmienna nBins.

clear all; close all; clc

%% Generacja sygnałów losowych

% Parametry próbkowania (niekonieczne przy sygnale losowym, ale dla porządku
% przyjmijmy 1kHz i czas obserwacji sygnału 5 sekund.)
fpr = 1000; dt = 1/fpr;     % częstotliwość próbkowania 

t0 = 5;  xM = 450;          % czas obserwacji sygnałów, alternatywnie t0 = 1 lub 0.1s. 

N = t0/dt;                  % liczba próbek
t = (0:N-1)*dt;             % wektor czasu

% typ = 1 - rozkład normalny, typ = 2 - rozkład Chi-kwadrat
typ = 1;

for h = 1 : 30
    
    switch(typ)
        case(1)
            u = 1; s = 0.5;         % u - wartość oczekiwana, s - odchylenie standardowe
            x = s*randn(1,N) + u;   % generujemy sygnał losowy o rozkładzie normalnym
            xMin = -2; xMax = 4; yMax = 450;    % ograniczenie osi (do rysunku)
            nBins = 50;
        case (2)
            L = 4;                % L - liczba stopni swobody
            x = chi2rnd(L, 1, N);   % generujemy sygnał losowy o rozkładzie chi-kwadrat
            xMin = -1; xMax = 10; yMax = 400;   % ograniczenie osi (do rysunku)
            nBins = 100; 
        otherwise
            return;
    end
    
    % Dla porządku obliczmy średnią i odchylenie standardowe
    uX=mean(x);
    sX=std(x);
       
    % Rysunki
    figure(1); subplot(211); plot(t,x);
    grid on; xlabel('Czas [s]'); ylabel('Wartość sygnału'); title(['Sygnał stochastyczny']);
    
    % Obliczanie histogramu (MATLAB histogram)
    subplot(212); histogram (x, nBins); title(['Histogram, średnia=' num2str(uX) '; odch. std.=' num2str(sX) ';']); 
    xlabel('Wartość sygnału'); ylabel('Liczność'); grid on;
    axis([xMin xMax 0 yMax]);   
    pause(1)
end

Rozkład normalny:

Rozkład Chi-kwadrat:

Co wynika z analizy rysunków? Najpierw wnioski ogólne:

• z histogramu możemy wnioskować jak często dana wartość występuje w sygnale, co jest podstawą określenia rodzaju rozkładu prawdopodobieństwa sygnału (o tym następnym razem) ,
• każda realizacja sygnału losowego ma inny kształt i inne wartości - sygnały w istocie są losowe,
• tym samym, kształt histogramu czyli częstości występowania poszczególnych wartości w sygnale, również się zmienia z realizacji na realizację (kolejno generowane sygnały), w sygnale o rozkładzie Chi-kwadrat nie występują liczby ujemne.
• kształty rozkładów wartości (histogramów) dla liczb normalnych i Chi-kwadrat różnią się, a w sygnale o rozkładzie Chi-kwadrat nie występują liczby ujemne,
• najczęściej występująca wartość dla rozkładu Normalnego to 1, a dla Chi-kwadrat to 2,
• w przypadku asymetrycznego rozkładu Chi-kwadrat wartość średnia (4) nie jest równa wartości najczęściej występującej (1).

Dużo cennych informacji!

Spróbujcie zmieniać czas obserwacji sygnału oraz liczbę przedziałów nBins i zaobserwujcie jak wpływa to na wynik obserwacji. Nie wdając się w analizę teoretyczną, zalecenia są takie, że przy analizę sygnałów losowych dąży się do obserwacji jak najdłuższego sygnału, a liczbę przedziałów nBins dobiera się tak, aby otrzymany obraz był czytelny.

Obiekt typu histogram

Wywołanie funkcji w następujący sposób:

h = historam(X, nBins);

powoduje utworzenie obiektu h typu histogram. Obiekt ten zawiera wszelkie informacje na temat obliczonego histogramu, które można łatwo odczytać za pomocą indeksowania z użyciem kropki, a ponadto daje to możliwości zmiany parametrów histogramu bez konieczności jego ponownego obliczania. Ilustruje to kolejny przykład:

clear all; close all; clc

%% Wczytanie danych i narysowanie histogramu

load carbig.mat     % Wczytujemy plik z danymi samochodów wyprodukowanych 
                    % w latach 1970 - 1982
                    
x = Horsepower;     % Interesuje nas histogram mocy pojazdów

nBins0=50;           % Liczba przedziałów histogramu
h = histogram(x,nBins0)  % h - obiekt typu histogram
xlabel('Moc silnika [KM]'); ylabel('Liczność'); grid on; 
disp(['Histogram z liczbą przedziałów ' num2str(nBins0)])
pause

%% Analiza histogramu i zmiana jego parametrów
disp('Liczność w poszczególnych przedziałach')
hVal = h.Values     % przez indeksowanie z kropką, mamy dostęp
                    % do poszczególnych wartości/parametrów
                    % histogramu. Values - to liczba zliczeń
                    % w poszczególnych przedziałach

disp('Liczba przedziałów:')
nBins = h.NumBins   % NumBins to liczba przedziałów
pause

% Liczbę przedziałów możemy zmieniać, już po wywołaniu funkcji 'histogram':

disp('Zwiększamy liczbe przedziałów: ')
h.NumBins = nBins + 10;
h.NumBins
pause

nBins = morebins(h)
pause

nBins = morebins(h)
pause

% Zamiast zmieniać liczbę przedziałów, można dostosować ich szerokość:
disp('Zamiast tego można zmieniać szerokość przedziału odwołując się tak: h.BinWidth')
wBin = h.BinWidth       % Obecna szerokość przedziału
h.BinWidth = 4 * wBin;  % Nowa szerokość przedziału (4 x stara)

Posługiwanie się obiektem histogram daje wiele możliwości i jest bardzo wygodne. Zamiast obliczać kolejne histogramy w celu otrzymania żądanego wyniku, można łatwo dostosowywać histogram do swoich potrzeb. Obiekt tego rodzaju, przechowuje również wyniki liczbowe dzięki temu można łatwo „wyciągnąć” dane liczbowe, na przykład liczby zliczonych wartości w danych przedziałach h.Values bądź zobaczyć jakie są granice przedziałów h.Binedges.

Funkcja histogram kryje w sobie dużo więcej niż z pozoru się wydaje. Można jej używać do wyznaczania funkcji gęstości prawdopodobieństwa sygnałów oraz analizy porównawczej różnych realizacji procesów losowych. Ale o tym już następnym razem. Tymczasem sami eksperymentujcie z funkcją MATLAB histogram.