% Napredne metode digitalne obradbe signala
% Kompresija audio signala

% Autor: Peter Skoda

% ************************************************************************
clear all; close all;

% *** Parametri **********************************************************

% Dekompozicija
wave_ime='sym4'; % MATLAB naziv wavelet funkcije

% Dodijela bitova i kvantizacija
b_sr=3; % prosjecni broj bitova po uzorku ulaznog signala

% *** Ucitavanje signala *************************************************
[s,fs,bit]=wavread('jamiroquai_01.wav');
s_l=length(s);

% *** Dekompozicija ******************************************************
% Radi se wavelet paket dekompozicija ortogonalnim waveletom. Rezultat
% dekompozicije je 11 pojasa, cvorovi dekompozicijskog stabla: [5 0],
% [5 1], [5 2], [5 3], [4 2], [4 3], [4 4], [4 5], [3 3], [2 2], [2 3].

% *** Dekompozicija *****
D=wpdec(s,2,wave_ime);
D=wpsplt(D,[2 0]);
D=wpsplt(D,[2 1]);
D=wpsplt(D,[3 0]);
D=wpsplt(D,[3 1]);
D=wpsplt(D,[3 2]);
D=wpsplt(D,[4 0]);
D=wpsplt(D,[4 1]);

% Koeficijenti dekompozicije, varijanca, broj koeficijenata
j=get(D,'tn');
for i=1:11
    d{i}=wpcoef(D,j(i)); % koeficijenti
    d_varijanca(i)=var(d{i}); % varijanca
    d_brkoef(i)=length(d{i}); % broj koeficijenata
end

% *** Dodijela bitova ****************************************************
% Svakom pojasu dodijeli se broj bitova prema pojednostavljenom "greedy"
% algoritmu. Podrazumjeva se da su normalizirane raspodijele vrijednosti
% koeficijenata jednake za svaki pojas.

broj_bitova=b_sr*length(s); % ukupni broj bitova na raspolaganju
% *** Dodijela *****
b=bit2(broj_bitova,d_varijanca,d_brkoef);

% *** Kvantizacija *******************************************************
% Vrsi se uniformna kvantizacija.

% *** Provjera bitova *****
% Ako je nekom pojasu dodijeljenu nula bitova, pojas se ne prenosi, tj. pri
% rekonstrukciji se uzima da su svi koeficijenti u tom pojasu jednaki nuli.
b_nz=find(b);

% *** Uniformna kvantizacija *****
for i=1:11,
    if b(i)==0;
        d_kv{i}=zeros(1,d_brkoef(i));
    else
        d_kv{i}=nckvantvar(d{i},b(i));
    end
end

% *** Rekonstrukcija *****************************************************

% *** Cijeli signal *****
A=get(D,'terNI');
A{2}=[d_kv{1};d_kv{2};d_kv{3};d_kv{4};d_kv{5};d_kv{6}; ...
        d_kv{7};d_kv{8};d_kv{9};d_kv{10};d_kv{11}].';
D=set(D,'terNI',A);
s_rec_sve=wprec(D);

% *** Sum pojedinih pojasa *****
A=get(D,'terNI');
A_ind=[0;cumsum(A{1}(:,1))];
for i=1:11
    A{2}=zeros(1,length(A{2}));
    A{2}(A_ind(i)+1:A_ind(i+1))=d_kv{i}-d{i};
    D=set(D,'terNI',A);
    s_rec_sum{i}=wprec(D);
end

% *** Aftermath **********************************************************

% Odnos signal/sum
s_raz=s_rec_sve-s;
SNR=10*log10(var(s)/( s_raz.'*s_raz/length(s) ));
clear s_raz;

% Srednja snaga pogreske po pojasu
for i=1:11
    N_rec(i)=s_rec_sum{i}.'*s_rec_sum{i}/length(s_rec_sum{i});
end

% *** Kraj ***************************************************************