Sygnały losowe i szumy
Autorzy: Przemysław Korohoda, Adam Borowicz, Krzysztof Duda
Pliki: rozdzial_03.zip
Przykład 3.1
Przykład (plik przyklad_3_1.m) zawiera kompletny plik skryptowy pokazujący jak samodzielnie eksperymentować z danymi pseudolosowymi na przykładzie rozkładu normalnego, a szczególnie jak przygotować obliczenia, aby następnie graficznie: a) porównać histogram z rozkładem teoretycznym, b) zweryfikować obciążenia estymatora odchylenia standardowego, c) zbadać rozkłady losowe estymatorów – wartości średniej i odchylenia standardowego.% Przykład 3.1:
% m-plik: przyklad_3_1.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 3. Sygnały losowe i szumy.
% Autorzy rozdziału: Przemysław Korohoda, Adam Borowicz, Krzysztof Duda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK.
clc; clear; close all; % posprzątanie "biurka";
N=2^10; K=2^14; % N-dlugość ciągu; K-liczba realizacji;
mu=3; sigma=4; % założone, czyli dokładne, parametry rozkładu;
X=randn(K,N)*sigma+mu; % generowanie kompletu danych (ciągi ergodyczne);
mx=mean(X); sx1=std(X); sx2=std(X,1); % szacowanie parametrów dla kolejnych "n";
[hx,bx]=hist(X(:),100); dx=bx(2)-bx(1); % wyliczanie histogramu i szerokości jednego "słupka";
[hmx,bmx]=hist(mx,100); % wyliczanie histogramu dla mx;
[hsx,bsx]=hist(sx1,100); % wyliczanie histogramu dla sx1 (estymator nieobciążony);
msx1=mean(sx1), msx2=mean(sx2), % wartości średnie dla rozkładów sx1 i sx2;
smx=std(mx), smx_ref=sigma/sqrt(K), % porównanie odchylenia standardowego dla mx z teorią;
figure(1); % przykład "surowego" opracowania graficznego:
subplot(311); bar(bx,hx/(sum(hx)*dx),1); grid on; hold on;
plot(bx,normpdf(bx,mu,sigma),'r-');
legend('rozklad teoretyczny','skalowany histogram');
subplot(312); bar(bmx,hmx,1); grid on;
legend('histogram wartosci sredniej');
subplot(313); bar(bsx,hsx,1); grid on;
legend('histogram odchyl. stand.');
% KONIEC PRZYKŁADU 3.1.
Przykład 3.2
Przykład (plik przyklad_3_2.m) zawiera kompletny plik skryptowy pokazujący jak samodzielnie eksperymentować z funkcją autokorelacji. Dane są generowane – zgodnie wybraną wersją – na 3 sposoby, a następnie przetwarzane zbiorowo. Na trójwymiarowych wykresach można zaobserwować: a) jakie są kształty i zmienność rozkładów, pokazanych w postaci histogramów, w kolejnych chwilach n, b) jakie są kształty i zmienność funkcji autokorelacji wyznaczanej dla kolejnych realizacji procesu, c) można porównać efekt zastosowania korelacji bez korekty obciążenia i z korektą (przesuwając znak komentarza na początku odpowiedniej linii do linii powyżej). Przykład zawiera możliwość wyboru jednego z trzech sposobów generowania danych pseudolosowych: 1) ergodyczny gaussowski proces losowy o zadanych parametrach, 2) proces dający realizacje kosinusoidalne o pseudolosowej fazie początkowej, 3) gaussowski proces losowy o skorelowanych wzajemnie kolejnych wartościach x[n].% Przykład 3.2:
% m-plik: przyklad_3_2.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 3. Sygnały losowe i szumy.
% Autorzy rozdziału: Przemysław Korohoda, Adam Borowicz, Krzysztof Duda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK.
clc; clear; close all; % Patrz Przykład 3.1;
N=2^8; K=2^10; % Patrz Przykład 3.1;
wersja=1; % wybór sposobu generowania danych (1,2 lub 3);
switch wersja,
case 1, mu=3; sigma=4; X=randn(K,N)*sigma+mu;
case 2, X=cos(ones(K,1)*2*pi*4*(0:N-1)/N + rand(K,1)*ones(1,N)*2*pi);
case 3, mu=3; sigma=4; X=randn(K,N)*sigma+mu;
M=5; h=rand(1,M); X=filter(h,1,X);
end;
[H,b]=hist(X,100); % wyznaczone zbiorczo histogramy (dla kolejnych realizacji);
opcja=1; % wybór sposobu wyliczania;
switch opcja,
case 1, for k=1:K, CX(k,:)=xcorr(X(k,:),X(k,:)); end; % wyliczamy rxx bez korekty;
case 2, m=-(N-1):N-1; % wariant z korektą;
for k=1:K, CX(k,:)=xcorr(X(k,:),X(k,:))./(N-abs(m)); end;
end;
figure(1);
mesh(0:N-1,b,H); title('wykres 3D z histogramami dla kolejnych realizacji procesu');
axis tight; xlabel('numer realizacji'); ylabel('x');
figure(2);
mesh(CX); title('wykres r_x_x dla kolejnych realizacji procesu');
axis tight; xlabel('m'); ylabel('numer realizacji');
% KONIEC PRZYKŁADU 3.2.
Przykład 3.3
W przykładzie (plik przyklad_3_3.m) pokazano jak (identyfikację odpowiednich elementów programu należy potraktować jako formę zadania): a) zweryfikować podobieństwa między splotem i funkcją korelacji, b) wyliczyć charakterystyki częstotliwościowe sygnału na rożne opisane dotąd sposoby. W tym przykładzie porównawcze studium otrzymanych wyników, w tym w formie graficznej, pozostawiono do samodzielnej realizacji. Należy zwrócić uwagę, że drugi wariant generowania danych dotyczy ciągów zespolonych.% Przykład 3.3:
% m-plik: przyklad_3_3.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 3. Sygnały losowe i szumy.
% Autorzy rozdziału: Przemysław Korohoda, Adam Borowicz, Krzysztof Duda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./PK.
clc; clear; close all;
K=2^10;
N=2^14; % Patrz przykład 3.1
F=(0:N-1)/N; % przygotowanie osi unormowanej częstotliwości;
wersja=1; % wybór wersji danych (wybór procesu losowego x);
switch wersja,
case 1, mu=0; sigma=4; x=randn(1,N)*sigma+mu; k=1:N/2+1; a=1;
case 2, muR=0; sigmaR=1; muI=0; sigmaI=1;
x=randn(1,N)*sigmaR+muR + j*(randn(1,N)*sigmaI+muI); % proces zespolony!!! Widmo nie będzie miało charakterystycznych symetrii;
k=1:N; a=1;
end;
y=rand(1,N); % realizacja procesu y;
c1=xcorr(y,x); c2=conv(y,conj(x(end:-1:1))); % wyznaczamy funkcję korelacji dwiema metodami;
disp('Roznica po wyznaczeniu funkcji korelacji dwiema metodami:');
err=max(abs(c1-c2)), % porównujemy wyniki;
X=fft(x); % DFT ciągu x;
[P1,W1]=periodogram(x,hamming(N),K); F1=W1/(2*pi); % periodogram z (pojedynczym) oknem Hamminga;
[P2,W2]=pwelch(x,K); F2=W2/(2*pi); % periodogram wyznaczony metodą Welcha;
figure(1);
plot(F1,P1,'ro-'); grid on; hold on; axis tight;
plot(F2,P2,'b.-');
xlabel('F'); ylabel('Widmowa gestosc mocy (PSD)'); legend('a) periodogram','b) pwelch');
title('Dwa periodogramy - konfrontacja wyników zastosowania dwoch funkcji');
% KONIEC PRZYKŁADU 3.3.
Przykład 3.4
Przykład 3.4 (plik przyklad_3_4.m) przedstawia kompletny plik skryptowy zawierający implementację generatora szumu kolorowego w pakiecie MATLAB, wykorzystując dolnoprzepustowy filtr typu IIR. Program umożliwia porównanie charakterystyk filtra z otrzymanym widmem PSD, do czego są konieczne pewne przeliczenia - warto się zastanowić z czego te przeliczenia wynikają.% Przykład 3.4:
% m-plik: przyklad_3_4.m
%
% Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014.
% Rozdział 3. Sygnały losowe i szumy.
% Autorzy rozdziału: Przemysław Korohoda, Adam Borowicz, Krzysztof Duda.
%
% Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl
%
% Opracowanie przykładu: kwiecień 2013 r./AB/PK.
clc; close all; clear;
N= 2^14; % dlugość przykładowej realizacji szumu;
fpr = 8000; % zakładamy częstotliwość próbkowania;
a = 0.9; % parametr filtra dolnoprzepustowego;
b1=[1-a]; a1=[1,-a]; % współczynniki filtra (patrz rozdział 4);
[H1,W1]=freqz(b1,a1,1000); % chki-filtra (patrz rozdział 4)
f1=fpr*W1/(2*pi);
sigma = 3; % odchylenie standard. dla PDF szumu białego;
x=randn(1,N)*sigma; % generujemy relizację procesu (szum biały);
y=filter(b1,a1,x); % modyfikujemy pasmo szumu przez filtrację;
[P,W2]=pwelch(y,256); % wyznaczamy periodogram;
f2=fpr*W2/(2*pi);
Pbis=P*pi/(sigma^2); % przeliczenia (pytanie-zadanie: dlaczego tak?);
figure(1);
plot(f2,10*log10(Pbis),'r.-'); grid on; hold on;
plot(f1,20*log10(abs(H1)),'c.-');
xlabel('f [Hz]'); ylabel('[dB]');
title('Porownanie ch-ki amplitudowej filtra "kolorujacego" i PSD otrzymanego szumu');
legend('PSD szumu kolorowego','ch-ka filtra (odpowiednio przeliczona)');
% KONIEC PRZYKŁADU 3.4.
Przykład 3.5
W przykładzie (plik przyklad_3_5.m) zawarto kompletny plik skryptowy weryfikujący zastosowania wzorów (3.45) i (3.46). Pod koniec zaproponowano sprawdzanie polegające na tym, że sygnał „wybielony blokami” został połączony w całość i podzielony na nowe bloki, tym razem na zakładkę – i ponownie wyliczono macierz kowariancji, aby ją porównać z zależnością (3.46).% Przykład 3.5: % m-plik: przyklad_3_5.m % % Cyfrowe przetwarzanie sygnałów. Podstawy, multimedia, transmisja. PWN, Warszawa, 2014. % Rozdział 3. Sygnały losowe i szumy. % Autorzy rozdziału: Przemysław Korohoda, Adam Borowicz, Krzysztof Duda. % % Więcej na stronie internetowej: http://teledsp.kt.agh.edu.pl % % Opracowanie przykładu: kwiecień 2013 r./PK. clc; clear; close all; N=2^10; M=8; K=N/M; % mamy zadbać, żeby się zgadzało (tj. K - liczba całkowita); sigma=3; mu=5; x=randn(1,N)*sigma+mu; % generujemy szum; x=filter([1,1,1]/3,1,x); % kolorujemy szum; X=reshape(x,M,K); % układamy bloki sygnału w macierzy X; m1=mean(X')'; X1=X-m1*ones(1,K); % odejmujemy średnie od współrzędnych; Rx= X1*X1', % wyznaczamy macierz kowariancji; Y=(Rx)^(-1/2)*X1; % wybielamy; Ry=Y*Y' , % sprawdzamy, czy wyszło; y=Y(:)'; % ... dokładniej sprawdzamy; for m=1:M, Y1(m,:)=y(m:end-M+m); end; % tworzymy bloki, ale na "zakładkę"; Ry_test=(1/M)*Y1*Y1', % ponowna weryfikacja; % KONIEC PRZYKŁADU 3.5.