Lab 5 — Τυχαία σήματα & θόρυβος

Στη θεωρία είδαμε τυχαίες μεταβλητές, τυχαίες διαδικασίες, PSD και λευκό θόρυβο. Εδώ τα φέρνουμε όλα στο MATLAB.

Γεννήτριες τυχαίων αριθμών

Παράδειγμα:

x = rand(100, 1);      % 100 τιμές U(0,1)
stem(x)                 % διακριτές τιμές
figure
hist(x, 100)           % ιστόγραμμα με 100 bins

Στατιστικά: mean, var, std

x = randn(1, 10000);            % 10000 τιμές N(0,1)
mean(x)                          % ≈ 0
var(x)                           % ≈ 1
std(x)                           % ≈ 1
median(x)                        % ≈ 0

Για κανονική N(μ, σ²) με μ=2, σ=4:

y = normrnd(2, 4, 1, 10000);
mean(y)     % ≈ 2
std(y)      % ≈ 4
var(y)      % ≈ 16

Σύγκριση δύο διανυσμάτων — covariance και correlation:

% Covariance
cov(x, y)
% Επιστρέφει 2x2 πίνακα: [var(x), cov(x,y); cov(x,y), var(y)]

% Correlation coefficient
corrcoef(x, y)
% Ομοίως 2x2 πίνακας, με ρ μεταξύ -1 και +1 στο off-diagonal

Τυχαία σήματα — generation patterns

Από τη θεωρία ξέρουμε ότι ένα τυχαίο σήμα έχει πραγματοποιήσεις (realizations). Στο MATLAB φτιάχνουμε πολλές πραγματοποιήσεις δοκιμάζοντας ξανά και ξανά.

Σήμα με τυχαίο πλάτος

NUM_REAL = 4;                            % πόσες πραγματοποιήσεις
N = 1024;                                 % δείγματα ανά πραγματοποίηση
t = 0:(1/N):(1 - 1/N);                   % ένα δευτερόλεπτο

figure;
for k = 1:NUM_REAL
    A = randi(4);                         % τυχαίο A από {1,2,3,4}
    x = A * cos(2*pi*4*t);
    subplot(NUM_REAL, 1, k);
    plot(t, x);
    title(sprintf('Πραγματοποίηση %d, A = %d', k, A));
end

Σήμα με τυχαία φάση

for k = 1:NUM_REAL
    theta = unifrnd(-pi, pi);             % τυχαία φάση
    x = cos(2*pi*4*t + theta);
    subplot(NUM_REAL, 1, k);
    plot(t, x);
end

Σήμα με τυχαία συχνότητα

for k = 1:NUM_REAL
    f = randi(4);                         % τυχαία ακέραια συχνότητα 1-4
    x = cos(2*pi*f*t);
    subplot(NUM_REAL, 1, k);
    plot(t, x);
end

Λευκός Gaussian θόρυβος (WGN)

Από τη θεωρία: ο λευκός θόρυβος έχει σταθερή PSD σε όλες τις συχνότητες (διπλάσιο νόημα: «λευκός» όπως το λευκό φως). Οι τιμές του είναι IID, και αν είναι Gaussian κάθε δείγμα είναι N(0, σ²).

N = 1024;
sigma = 1;
n = sigma * randn(1, N);                 % WGN με τυπική απόκλιση σ

figure;
subplot(2,1,1);
plot(n)
title('White Gaussian Noise')

subplot(2,1,2);
hist(n, 50)
title('Ιστόγραμμα — πρέπει να μοιάζει με καμπύλη Gauss')

Σήμα + θόρυβος

fs = 1000;                                % συχνότητα δειγματοληψίας
t = 0:1/fs:0.5;                           % 0.5 s
A1 = 3; f1 = 270;
A2 = 8; f2 = 175;
signal = A1*cos(2*pi*f1*t) + A2*cos(2*pi*f2*t);

% Προσθήκη θορύβου με μέση τιμή 3 και διασπορά 25
noise = 5*randn(size(t)) + 3;             % σ=5 ⟹ σ²=25, μ=3
X = signal + noise;

plot(t, X)
title('Σήμα + θόρυβος — φαίνεται "βρώμικο"')

Το ωραίο: στο γράφημα φαίνεται απλώς θόρυβος. Στο πεδίο των συχνοτήτων όμως, οι 2 τόνοι (στις 175 και 270 Hz) εμφανίζονται καθαρά. Αυτό μας το δείχνει το PSD.

Εκτίμηση PSD με periodogram

Από τη θεωρία: PSD = μετασχηματισμός Fourier της αυτοσυσχέτισης. Στο MATLAB υπάρχει έτοιμη εκτιμήτρια — το περιοδόγραμμα.

[S, w] = periodogram(X, [], 1024, 1);     % () placeholder, 1024 FFT points, fs=1
plot(w, S)
xlabel('\omega'), ylabel('PSD'), grid on
title('Εκτίμηση PSD')

Παράδειγμα — δύο ημίτονα μέσα σε θόρυβο

fs = 1000;
t = 0:1/fs:0.5;
X = 5*randn(size(t)) + 3 + 3*cos(2*pi*270*t) + 8*cos(2*pi*175*t);

[S, w] = periodogram(X, [], 1024, 1);

subplot(2,1,1)
plot(t, X), title('Σήμα στο χρόνο — δεν διακρίνονται οι τόνοι')

subplot(2,1,2)
plot(-2*w, S, 2*w, S), grid on
title('PSD — οι τόνοι ξεχωρίζουν καθαρά στις 175 και 270 Hz')
xlabel('\omega/\pi')

Θόρυβος μέσα από φίλτρα

Όταν φιλτράρουμε λευκό θόρυβο, διαμορφώνουμε το φάσμα του θορύβου. Λευκός θόρυβος = επίπεδη PSD. Μετά από LP φίλτρο → η PSD γίνεται «χαμηλοπερατή» (μόνο χαμηλές συχνότητες απομένουν).

% Δημιουργία λευκού θορύβου
fs = 1000;
t = 0:1/fs:0.5;
X = 2*randn(size(t)) + 4;                 % μ=4, σ=2

% Κατασκευή χαμηλοπερατού φίλτρου τάξης 150, συχνότητα αποκοπής Wn=0.5
h = fir1(150, 0.5);                       % FIR filter
Y = filter(h, 1, X);                      % πέρασμα μέσα από φίλτρο

% PSD πριν και μετά
[SX, w] = periodogram(X, [], 1024, 1);
[SY, w] = periodogram(Y, [], 1024, 1);

subplot(2,1,1)
plot(-2*w, SX, 2*w, SX), grid on
title('PSD εισόδου — επίπεδη (λευκός θόρυβος)')

subplot(2,1,2)
plot(-2*w, SY, 2*w, SY), grid on
title('PSD εξόδου — αποκόπηκαν οι υψηλές συχνότητες')
xlabel('\omega/\pi')

Για υψηπερατό φίλτρο, αρκεί:

h = fir1(150, 0.5, 'high');

Συνάρτηση αυτοσυσχέτισης

Από τη θεωρία: η αυτοσυσχέτιση R_X(τ) = E[X(t)X(t+τ)] είναι το βασικό εργαλείο για στατιστικές διεργασίες. Στο MATLAB υπολογίζεται με xcorr.

X = randn(1, 10000);                      % WGN
r = xcorr(X, 'biased');                   % unbiased option = κανονικοποιημένη

% Εστίαση γύρω από το lag = 0
N = 100;
mid = length(r)/2;                        % κέντρο
stem(-N:N, r(mid-N:mid+N))
title('Αυτοσυσχέτιση WGN — μοιάζει με δέλτα(τ)')
xlabel('lag τ'), ylabel('R(\tau)')

Για λευκό θόρυβο, η αυτοσυσχέτιση είναι σχεδόν 0 παντού εκτός του τ = 0, όπου έχει αιχμή ίση με την διασπορά. Αυτό ΕΙΝΑΙ η μαθηματική έκφραση του «άσχετα μεταξύ τους δείγματα».

Ασκήσεις δοκιμής

Quick recap

• Γεννήτορες: rand (ομοιόμορφη), randn (κανονική), randi (διακριτή ομοιόμορφη), normrnd, unifrnd για πιο γενικά.
• Στατιστικά: mean, var, std, cov, corrcoef.
• WGN: sigma * randn(1, N) δίνει λευκό Gaussian θόρυβο με τυπική απόκλιση σ.
• PSD: periodogram(X, [], NFFT, fs) εκτιμά το φάσμα.
• Φιλτράρισμα θορύβου: fir1(N, Wn) κατασκευάζει το φίλτρο, filter(h, 1, X) φιλτράρει.
• Αυτοσυσχέτιση: xcorr(X, 'biased').