TeleDSP

Estymacja i redukcja zakłóceń w sygnale mowy

Autorzy: Adam Borowicz

Pliki: rozdzial_13.zip

Przykład 13.1

Przykład (plik przyklad_13_1.m) realizuje uśrednianie wykładnicze widma mocy w ramkach . Argumentami wejściowymi funkcji są: wektor zawierający próbki sygnału, parametr uśredniający oraz indeks próbki częstotliwościowej, dla której są wyświetlane rezultaty. Zasadnicze znaczenie ma tu parametr uśredniający alfa, od którego zależy nie tylko wariancja estymatora, ale również jego obciążenie.

function expavg(x, alpha, k)

N = 256; 
K = 128;
x = x(:); 
l = 1; 
Pdd(:,l) = zeros(N,1);

for n = 1:N-K:length(x)-N,
    l=l+1;
    % krótkookresowe widmo mocy
    X2(:,l) = abs(fft(x(n:n+N-1))).^2;
    % uśrednianie wykładnicze widma mocy w ramkach
    Pdd(:,l) = alpha * Pdd(:,l-1) + (1-alpha) * X2(:,l);
end;

plot(10*log10( X2(k,:)+1e-6), ’b:’), hold on;
plot(10*log10(Pdd(k,:)+1e-6), ’r-’);

Przykład 13.2

Przykład (plik przyklad_13_2.m) przedstawia implementację detektora energetycznego. Ze względu na nieznany początkowy poziom energii szumu algorytm może generować błędy, dlatego też do jego poprawnej pracy jest wymagana początkowa faza adaptacji, umożliwiająca dostrojenie parametrów początkowych do aktualnej charakterystyki zakłóceń.

function ym = edet(x, alpha, lambda)

N = 256;   % długość ramki
K = 128;   % nakładanie się ramek
adpt = 40; % czas adaptacji (liczba ramek)
x = x(:);
m = 0;

for n = 1:N-K:length(x)-N,
    m = m + 1;

    xw = x(n:n+N-1);                % obliczanie energii w m-tej ramce
    ex(m) = xw'*xw;

    if m > adpt & ex(m) > thr(m-1), % wypowiedź
        ym(m) = 1;
        thr(m) = thr(m-1);
    else                            % przerwa w wypowiedzi
        ym(m) = 0;
        if m == 1, ed = ex(m); ed2 = ed^2; end;
        ed  = alpha*ed  + (1-alpha)*ex(m);           % est. wartości oczek.
        ed2 = alpha*ed2 + (1-alpha)*ex(m)^2;         % est. wariancji
        thr(m) = ed + lambda*sqrt(abs(ed2 - ed2^2)); % próg decyzyjny
    end;
end;

Przykład 13.3

Przykład (plik przyklad_13_3.m) zawiera uproszczony fragment kodu, realizujący mechanizm korekcji błędów przypadkowych.

count1 = 0;
count2 = 0;
M = 3;

for m = 1:liczba_ramek,
    I = ...    % obliczanie decyzji detektora (1 - mowa, 0 - przerwa)
    if I == 1, 
        % korekcja błędów w trakcie wypowiedzi
        count2 = 0; count1 = count1 + 1;
        if count1 > M, ym(m) = 1; else ym(m) = ym(m-1); end;
    else       
        % korekcja w przerwach
        count1 = 0; count2 = count2 + 1;
        if count2 > M, ym(m) = 0; else ym(m) = ym(m-1); end;
    end;
end;

Przykład 13.4

Przykład (plik przyklad_13_4.m) przedstawia implementację estymatora MCRA w dla l-tej ramki sygnału. Funkcja wymaga przekazania estymaty widmowej mocy sygnału zakłóconego w postaci wektora Pxx. Konieczna jest również inicjalizacja struktury parameters, zawierającej oszacowania parametrów z poprzedniej ramki.

function parameters = MCRA_estimation(Pxx, parameters)

alpha = 0.95;
beta = 0.1;
delta = 2.5;
L = 50;

% oszacowania z poprzedniej ramki
pkl = parameters.pkl;   % szacunkowe prawdopodob. wypowiedzi
l = parameters.l;       % indeks ramki
Pmin = parameters.Pmin; % estymata minimum statystyki
Ptmp = parameters.Ptmp; % wektor pomocniczy

if rem(l,L)==0
    Pmin = min(Ptmp,Pxx);
    Ptmp = Pxx;
else
    Pmin = min(Pmin,Pxx); Ptmp = min(Ptmp,Pxx);
end

Srk = Pxx./Pmin;
Ikl = zeros(length(Srk),1);
Ikl(find(Srk > delta)) = 1;
pkl = beta*pkl + (1-beta) * Ikl;       % prawdopodob. wypowiedzi
alpha_t = alpha + (1-alpha)*pkl;       % parametr uśredniający
Pdd = alpha_t.*Pdd + (1-alpha_t).*Pxx; % estymacja PSD szumu

% zapis parametrów dla ramki l+1
parameters.pkl = pkl; 
parameters.l = l+1;
parameters.Pdd = Pdd;
parameters.Pmin = Pmin;
parameters.Ptmp = Ptmp;

Przykład 13.5

Program w języku MATLAB (plik przyklad_13_5.m), realizujący metodę różnic widmowych. W algorytmie tym przyjęto dla uproszczenia, że początkowe okna sygnału zakłóconego zawierają wyłącznie szum. Widmo szumu jest estymowane dla pierwszych 10 ramek sygnału, z wykorzystaniem średniej arytmetycznej.

function y = ss_nr(x)

x = x(:);     % wymuszenie notacji kolumnowej
N = 256;      % długość ramki (w próbkach), musi być parzysta!
K = 128;      % nakładanie się ramek (l. próbek)
W = hann(nw); % funkcja okna analizy
a = 1;        % typ metody: 1 - MSS, 2 - PSS
lmbd = 3;     % współczynnik kompensacji
N2 = N/2+1; 
y = zeros(size(x)); 
Da = zeros(N2,1); 
Xa = zeros(N2,1);

for n = 1:N-K:length(x)-N,
    % analiza
    sel = n:n+N-1;
    X = fft(x(sel).*W);
    Xa = abs(X(1:N2)).^a;
    % estymacja szumu w oparciu o 10 początkowych ramek
    if n < 10, Da = ((n-1)*Da + Xa)/n; end;
    % ważenie widmowe
    Gss = max(1 - lmbd*(Da./(Xa+eps)), 0).^(1/a);
    Y = X(1:N2) .* Gss;
    % synteza
    Y(N2+1:N) = flipud(conj(Y(2:N2-1)));   % odtworzenie symetrii DFT
    y(sel) = y(sel)+real(ifft(Y))./max(W); % nakładanie z dodawaniem
end;

subplot(2,1,1); 
    plot(x);
subplot(2,1,2); 
    plot(y);

Przykład 13.6

W przykładzie (plik przyklad_13_6.m) przedstawiono implementację estymatora TDC w języku MATLAB. Dla uproszczenia przyjęto założenie, że sygnałem zakłócającym jest szum biały, jednak w ogólnym przypadku należałoby przeprowadzić wstępne wybielanie sygnału zakłóconego. Ponadto w celu przyśpieszenia obliczeń, sąsiadujące ze sobą wektory x(l) są grupowane w kolumnach macierzy trajektorii, po czym przemnażane przez tę samą macierz KLT. Uzyskane współczynniki są modyfikowane przy użyciu diagonalnej macierzy wagowej i transformowane z powrotem do dziedziny czasu. Ramki sygnału uzdatnionego są uzyskiwane z wynikowej macierzy trajektorii poprzez uśrednianie diagonalne.

function y = tdc_nr(x)

N = 32;          % wymiar przestrzeni sygnału
M = 256;         % długość ramki (l. próbek)
K = 128;         % nakładanie się ramek (l. próbek)
Lmbd_thr = 3e-3; % parametr progowy dla wartości własnych
mu = 5;          % mnożnik Lagrange'a
x = x(:);
y = zeros(size(x));
W = bartlett(M);

% estymacja wariancji szumu na podstawie 1000 pierwszych próbek sygnału x
sigma2_d = var(x(1:1000));

for n = 1:M-K:length(x)-N, % przetwarzanie w ramkach
    
    sel = n:n+M-1;
    Xt = buffer(x(sel),N,N-1,’nodelay’); % mac. trajektorii
    Cx = (Xt*Xt’)./size(Xt,2);           % empiryczna mac. kowariancji
    Cy = Cx - sigma2_d*eye(N);           % est. mac. kowariancji mowy
    [Uy, Lmbd_y] = eig(Cy);              % est. mac. KLT i war. własnych
    Lmbd_y = diag(Lmbd_y);
    I = find(Lmbd_y > Lmbd_thr);         % progowanie war. własnych
    
    g_tdc = zeros(N,1);
    g_tdc(I) = Lmbd_y(I)./(Lmbd_y(I) + mu*sigma2_d); % wsp. filtru TDC
    Yt = Uy * (diag(g_tdc) * (Uy’ * Xt));            % filtracja
    
    % uśrednianie diagonalne (synteza sygnału)
    Y = zeros(M,1);
    w = zeros(M,1); 
    i = 1:N;
    for m = 1:size(Yt,2),
        Y(i) = Y(i) + Yt(:,m);
        w(i) = w(i) + 1;
        i = i+1;
    end;
    
    % nakładanie (uśrednianie) ramek z dodawaniem
    y(sel) = y(sel) + (Y./w).*W; 
end;

subplot(2,1,1); 
    plot(x);
subplot(2,1,2); 
    plot(y);