TeleDSP

Kompresja sygnałów fonicznych

Autorzy: Maciej Bartkowiak

Pliki: rozdzial_15.zip

Przykład 15.1

Jest to uproszczona implementacja modelu kodeka (plik przyklad_15_1.m), która dokonuje dekompozycji sygnału na wskazaną liczbę podpasm. Wartości próbek każdego podpasma są skalowane i kwantowane z określoną rozdzielczością bitową. Liczba bitów potrzebnych do zakodowania skwantowanych próbek naturalnym kodem binarnym jest szacowana na podstawie zakresu dynamicznego w każdym paśmie. Następnie sygnał jest rekonstruowany przez odwrócenie procesu. Program kodowanie_podpasmowe znajduje się w pliku plik kodowanie_podpasmowe.m.

% przykład 15.1: Demonstracja kodowania podpasmowego

clc; clear; close all
[x,Fs] = wavread('harpsichord.wav'); 	% Wczytanie sygnału z pliku wav
soundsc(x,Fs)                           % odsłuchanie
figure(1)
specgram(x,2048,Fs,2000)                % demonstracja spektrogramu
title('Sygnał oryginalny')

% Zakodowanie koderem podpasmowym z użyciem 8 podpasm i 6 bitów na próbkę
% UWAGA: przed skwantowaniem każde pasmo jest normalizowane do zakresu –1..1
[y2,bps] = kodowanie_podpasmowe(x,8,6);
soundsc(y2,Fs)
figure(2)
specgram(y2,2048,Fs,2000)
title(sprintf('Liczba bitów na próbkę: %1.2f\n',bps))

% Zakodowanie koderem podpasmowym z użyciem 32 podpasm i 8 bitów na próbkę
[y3,bps] = kodowanie_podpasmowe(x,32,5);
soundsc(y3,Fs)
figure(3)
specgram(y3,2048,Fs,2000)
title(sprintf('Liczba bitów na próbkę: %1.2f\n',bps))

% Zakodowanie koderem podpasmowym z użyciem 16 podpasm i zmiennej liczby bitów na próbkę
[y4,bps] = kodowanie_podpasmowe(x,16,[6 5 4 4 3 3 3 3 3 3 3 3]);
soundsc(y4,Fs)
figure(4)
specgram(y4,2048,Fs,2000)
title(sprintf('Liczba bitów na próbkę: %1.2f\n',bps))

Przykład 15.2

Szybka implementacja przekształcenia DCT-4 w formie mnożenia macierzy. Kod przykładu jest zapisany w pliku dct4.m.

function y = dct4(x)
% Szybkie obliczenie transformaty DCT-4 z bloku próbek
% y = dct4(x)
% x - pionowy wektor próbek
% y - pionowy wektor współczynników

x = x(:);
N = length(x);
if mod(N, 2)
    error('Rząd transformacji musi być parzysty');
end
N  = length(x);
[n,k] = meshgrid(0:(N-1),0:(N/2-1));
C = cos(pi*(2*n+1+N/2).*(2*k+1)/(2*N));
y = C*x;

Przykład 15.3

Szybka implementacja odwrotnego przekształcenia DCT-4 w formie mnożenia macierzy. Kod przykładu jest zapisany w pliku idct4.m.

function y = idct4(x)
% Szybkie obliczenie odwrotnej transformaty DCT-4 z bloku współczynników
%
% y = idct4(x)
% x - pionowy wektor współczynników 
% y - pionowy wektor próbek

x = x(:);
N = 2*length(x);
[k,n] = meshgrid(0:(N/2-1),0:(N-1));
C = cos(pi*(2*n+1+N/2).*(2*k+1)/(2*N));
y = 4*C*x/N;

Przykład 15.4

Demonstracja (plik przyklad_15_4.m), jak wygląda zrekonstruowany blok próbek sygnału po wykonaniu DCT-4 i odwrotnej DCT-4. Na wspólnym wykresie wyświetlimy sygnał oryginalny (niebieski), składniki aliasowe (zielone) oraz ich złożenie, czyli sygnał zrekonstruowany w 1 bloku (na czerwono).

% przykład 15.4

% Przygotowanie wektora kolumnowego próbek funkcji liniowo rosnącej
x = (1:100)';
% Obliczenie transformaty i odwrotnej transformaty
F = dct4(x);
y = idct4(F);

% Wykresy 
figure(1)
plot([x [-x(end/2:-1:1); x(end:-1:end/2+1)] y])
legend('Oryginalny','Nakładki','Wynik transformaty odwrotnej')

disp('Proszę nacisnąć dowolny klawisz');
pause

% Przygotowanie wektora próbek przykładowego sygnału
t = (1:1e3)';
x = exp(-3*t/1e3).*cos(2*pi*t.^2/1e5);
% Obliczenie transformaty i odwrotnej transformaty
F = dct4(x);
y = idct4(F);

% Wykresy 
plot([x [-x(end/2:-1:1); x(end:-1:end/2+1)] y])
legend('Oryginalny','Nakładki','Wynik transformaty odwrotnej')


disp('Proszę nacisnąć dowolny klawisz');
pause

% Przygotowanie wektora próbek przykładowego sygnału
t = linspace(-5,5,1e4)';
x = sinc(t);
% Obliczenie transformaty i odwrotnej transformaty
F = dct4(x);
y = idct4(F);

% Wykresy 
plot([x [-x(end/2:-1:1); x(end:-1:end/2+1)] y])
legend('Oryginalny','Nakładki','Wynik transformaty odwrotnej')

Przykład 15.5

Uproszczona implementacja kodera transformatowego, który dokonuje podziału sygnału x na bloki o długości N, z zakładką N/2 i oblicza transformatę MDCT. Współczynniki transformaty Fk są skalowane i kwantowane zgodnie z podanymi współczynnikami skalującymi Q, którymi można regulować stopień kompresji (jeden wspólny współczynnik lub wektor indywidualnych współczynników osobno dla każdego indeksu częstotliwości). Wyjściowa tablica symboli sym to całkowitoliczbowa reprezentacja skwantowanych współczynników transformaty. Wypadkowa wielkość zakodowanego strumienia (średnia liczba bitów na próbkę, bps) jest szacowana na podstawie zakresu dynamicznego symboli (osobno dla każdego indeksu częstotliwości). Implementacja kodera jest zapisana w pliku kodtr.m.

function [sym,bps] = kodtr(x,N,Q)
% Koder transformatowy
%
% [sym,bps] = kodtr(x, N, Q)
% x – sygnał wejściowy
%
% N – długość bloku próbek dla transformaty MDCT
% Q – współczynniki skalujące (jeden wspólny lub wektor indywidualnych współczynników)
% sym – tablica zakodowanych symboli
% bps – średnia liczba bitów zakodowanych danych na próbkę sygnału wejściowego

H = N/2;                        % przesunięcia kolejnych bloków
M = floor((length(x)-H)/H);      % liczba bloków
sym = zeros(H,M);               % tablica symboli danych zakodowanych
win = sin(pi*((0:(N-1))+0.5)/N)'; % okienko do transformaty MDCT

h_wbar = waitbar(0,'Przetwarzanie ramek', 'Name', 'Kodowanie transformatowe');
for m = 0:M-1
    waitbar(m/M,h_wbar);
    n0 = m*H + 1;               % początek bloku
    x0 = x(n0:n0+N-1);          % pobieranie bloku próbek
    x0 = x0.*win;               % okienkowanie
    Fk = dct4(x0);              % obliczenie transformaty
    Fkq = fix(Fk.*Q);           % kwantowanie współczynników
    sym(:,m+1) = Fkq;           % zapisanie do tablicy symboli
end
close(h_wbar);

% Oszacowanie wielkości strumienia danych
zakr = (max(sym')-min(sym'))';  % zakresy zmienności symboli
koszt = max(0,ceil(log2(zakr))); % szacunkowy koszt zakodowania symboli
bps = mean(koszt);              % średnia liczba bitów na próbkę

Przykład 15.6

Implementacja dekodera transformatowego zgodna z koderem z przykładu 15.5. Dekoder rekonstruuje sygnał y poprzez odwrócenie skalowania zastosowanego do współczynników transformaty przed kwantyzacją, zastosowanie odwrotnego przekształcenia MDCT i zsumowania odtworzonych bloków próbek z zakładką 50%. Implementacja funkcji dekodera znajduje się w pliku dektr.m.

function y = dektr(sym,N,Q)
% Dekoder transformatowy
% y = dektr(sym, N, Q)
%
% sym – tablica zakodowanych symboli
% N – długość bloku próbek dla transformaty MDCT
% Q – współczynniki skalujące (jeden wspólny lub wektor indywidualnych współczynników)
% y – sygnał odtworzony

H = N/2;                        % przesunięcia kolejnych bloków
M = size(sym,2);                % liczba bloków
L = H * (M+1);                  % szacowana długość sygnału zdekodowanego
y = zeros(L,1);                 % miejsce na próbki sygnału
win = sin(pi*((0:(N-1))+0.5)/N)'; % okienko do transformaty MDCT

h_wbar = waitbar(0,'Dekodowanie ramek', 'Name', 'Kodowanie transformatowe');
for m = 0:M-1
    waitbar(m/M,h_wbar);
    Fkq = sym(:,m+1);           % odczytanie kolejnego wektora symboli
    Fkr = Fkq./Q;               % odtworzenie współczynników transformaty
    y0 = idct4(Fkr);            % obliczenie odwrotnej transformaty
    y0 = y0.*win;               % ponowne okienkowanie
    n0 = m*H + 1;               % początek bloku    
    y(n0:n0+N-1) = y(n0:n0+N-1)+y0; % składanie bloków na zakładkę
end
close(h_wbar);

Przykład 15.7

Uproszczona implementacja bloku operacji TNS w koderze. Operacja jest powtarzana dla pojedynczych bloków współczynników transformaty kosinusowej, które odpowiadają blokom próbek, w których wykryto impuls (transient). Zmodyfikowane współczynniki zastępują oryginalne współczynniki w tym bloku. Kod funkcji znajduje się w pliku tns_pre.m.

function [Fk,R] = tns_pre(Fk,Pmax,Gmin)
% Fk = tns_pre(Fk,Pmax,Gmin)
%
% Fk - wektor współczynników transformaty MDCT
% Pmax - maksymalny rząd predykcji
% Gmin - minimalny zysk predykcji
% R - współczynniki PARCOR predyktora

Rmin = 0.05;

% wyznaczamy współczynniki predyktora
a = lpc(Fk,Pmax);

% przekształcamy do postaci PARCOR
R = poly2rc(a); 

% znajdujemy małe wartości współczynników
z = find(abs(R)>=Rmin,1,'last'); 
if ~isempty(z)
    % odrzucamy końcowe małe współczynniki zmniejszając rząd predyktora
    R(z+1:end)=[]; 
    a = rc2poly(R); % przeliczenie z PARCOR na LPC
end

% filtr analizy ("wybielanie" w czasie, 
% czyli wyrównywanie amplitudy sygnału)
X = filter(a,1,Fk); 

% obliczanie zysku predykcji
Gp = rms(Fk)/rms(X); 
if Gp >= Gmin % warunek aktywności narzędzia TNS
    % zastąpienie wybranych współczynników MDCT wartościami błędu predykcji
    Fk = X;     
else
    R = [];
end

Przykład 15.8

Uproszczona implementacja bloku operacji TNS w dekoderze. Operacja jest wykonywana dla tych samych pojedynczych bloków współczynników transformaty kosinusowej, które zostały zmodyfikowane po stronie kodera (informacja o włączeniu TNS wynika z niezerowego bloku współczynników predyktora powiązanego z danym blokiem transformaty). Rezultatem operacji jest odtworzenie właściwych współczynników, które reprezentują zakodowany sygnał. Kod funkcji jest zapisany w pliku tns_post.m.

function Fk = tns_post(Fk,R)
% Fk = tns_post(Fk,A)
%
% Fk - wektor współczynników transformaty MDCT
% R - współczynniki predyktora

% przeliczenie współczynników do postaci LP
a = rc2poly(R); 

% filtr syntezy (odtwarzanie przebiegu)
Fk = filter(1,a,Fk);