

% brisanje svih varijabli u workspaceu
clear
clc

% propuštanje zašumljenog signala kroz sustav

% generiranje testnog signala
% s ... originalni signal
% x ... zašumljeni signal
% n ... šum
s = wnoise(4, 10, 7);


% ucitavanje wav datoteke sa zvukom gitare
% [s, fs, nbits] = wavread('gitara.wav');
% uzimamo samo jedan kanal
% s = (s(1:end,1))';

% generiranje bijelog šuma koji dodajemo signalu
mu = 0;
sigma = 2;
duljina = length(s);
n = normrnd(mu, sigma, 1, duljina);

% zašumljeni signal na ulazu sustava
x = s + n;

% prikaz generiranih signala
figure('Name','Prikaz signala na ulazu')
subplot(3,1,1), plot(s);
title('s'); xlabel('uzorak'); ylabel('amplituda');
subplot(3,1,2), plot(n);
title('n'); xlabel('uzorak'); ylabel('amplituda');
subplot(3,1,3), plot(x);
title('x = s + n'); xlabel('uzorak'); ylabel('amplituda');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prva wavelet transformacija ulaznog signala i obrada sa hard thresholdom
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% zadajemo wavelet za prvu transformaciju
W_1 = 'sym5';

% dekompozicija u 5 razina
[y_1_C, y_1_L] = wavedec(x, 5, W_1);
y_1_a5 = appcoef(y_1_C, y_1_L, W_1, 5);
y_1_d5 = detcoef(y_1_C, y_1_L, 5);
y_1_d4 = detcoef(y_1_C, y_1_L, 4);
y_1_d3 = detcoef(y_1_C, y_1_L, 3);
y_1_d2 = detcoef(y_1_C, y_1_L, 2);
y_1_d1 = detcoef(y_1_C, y_1_L, 1);

% prikaz dekompozicije y_1
figure('Name','Prikaz dekompozicije y_1')
subplot(6,1,1); plot(y_1_a5); title('y_1 a5');
subplot(6,1,2); plot(y_1_d5); title('y_1 d5');
subplot(6,1,3); plot(y_1_d4); title('y_1 d4');
subplot(6,1,4); plot(y_1_d3); title('y_1 d3');
subplot(6,1,5); plot(y_1_d2); title('y_1 d2');
subplot(6,1,6); plot(y_1_d1); title('y_1 d1');

% procjena sigme iz medijana
% procjenjujemo je kao: 1.4826*medijan_apsolutne_vrijednosti
% ( vidi funkciju wnoisest )

K1 = 1.4826;
sigma_est = K1*median(abs(y_1_d1));

% IDEALNA PROCJENA SIGME ZA OVAJ SIGNAL JE IZ KOEFICIJENATA cD1, TJ.
% DETALJA U PRVOJ RAZINI RAZLAGANJA JER U TOM POJASU IMAMO SAMO ŠUM

% granica za threshold
% K2 odabiremo izmeðu 2 i 5
K2 = 4;
G = K2*sigma_est;

% obrada u domeni transformacije
% hard threshold sa granicom G
theta_1_a5 = wthresh(y_1_a5,'h',G);
theta_1_d5 = wthresh(y_1_d5,'h',G);
theta_1_d4 = wthresh(y_1_d4,'h',G);
theta_1_d3 = wthresh(y_1_d3,'h',G);
theta_1_d2 = wthresh(y_1_d2,'h',G);
theta_1_d1 = wthresh(y_1_d1,'h',G);


% prikaz koeficijenata nakon hard threshold-a
figure('Name','Prikaz koeficijenata theta_1')
subplot(6,1,1); plot(theta_1_a5); title('theta_1 a5');
subplot(6,1,2); plot(theta_1_d5); title('theta_1 d5');
subplot(6,1,3); plot(theta_1_d4); title('theta_1 d4');
subplot(6,1,4); plot(theta_1_d3); title('theta_1 d3');
subplot(6,1,5); plot(theta_1_d2); title('theta_1 d2');
subplot(6,1,6); plot(theta_1_d1); title('theta_1 d1');

theta_1_C = [theta_1_a5   theta_1_d5   theta_1_d4   theta_1_d3   theta_1_d2   theta_1_d1];
theta_1_L = y_1_L;

% inverzna transformacija
s_1 = waverec(theta_1_C, theta_1_L, W_1);

% prikaz rekonstruiranog signala nakon hard thresholda
figure('Name','Prikaz signala s_1')
subplot(4,1,1), plot(s);
title('s'); xlabel('uzorak'); ylabel('amplituda');
subplot(4,1,2), plot(x);
title('x = s + n'); xlabel('uzorak'); ylabel('amplituda');
subplot(4,1,3), plot(s_1);
title('s_1'); xlabel('uzorak'); ylabel('amplituda');
subplot(4,1,4), plot(s_1-s);
title('s - s1'); xlabel('uzorak'); ylabel('amplituda');



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% druga wavelet transformacija ulaznog signala i obrada sa empirijski
% dobivenim Wienerovim filtrom
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% druga wavelet transformacija ulaznog signala i izlaza iz prve wavelet
% transformacije

% zadajemo wavelet za drugu wavelet transformaciju W2
W_2 = 'sym4';

% dekompozicija u 5 razina wavelet transformacijom W2

% ulazni signal
[y_2_C, y_2_L] = wavedec(x, 5, W_2);
y_2_a5 = appcoef(y_2_C, y_2_L, W_2, 5);
y_2_d5 = detcoef(y_2_C, y_2_L, 5);
y_2_d4 = detcoef(y_2_C, y_2_L, 4);
y_2_d3 = detcoef(y_2_C, y_2_L, 3);
y_2_d2 = detcoef(y_2_C, y_2_L, 2);
y_2_d1 = detcoef(y_2_C, y_2_L, 1);

% aproksimacija ulaznog signala dobivena hard threshold-om
[theta_21_C, theta_21_L] = wavedec(s_1, 5, W_2);
theta_21_a5 = appcoef(theta_21_C, theta_21_L, W_2, 5);
theta_21_d5 = detcoef(theta_21_C, theta_21_L, 5);
theta_21_d4 = detcoef(theta_21_C, theta_21_L, 4);
theta_21_d3 = detcoef(theta_21_C, theta_21_L, 3);
theta_21_d2 = detcoef(theta_21_C, theta_21_L, 2);
theta_21_d1 = detcoef(theta_21_C, theta_21_L, 1);

% prikaz dekompozicije
figure('Name','Prikaz dekompozicije y_2 i theta_21')
subplot(6,2,1); plot(y_2_a5); title('y_2 a5');
subplot(6,2,3); plot(y_2_d5); title('y_2 d5');
subplot(6,2,5); plot(y_2_d4); title('y_2 d4');
subplot(6,2,7); plot(y_2_d3); title('y_2 d3');
subplot(6,2,9); plot(y_2_d2); title('y_2 d2');
subplot(6,2,11); plot(y_2_d1); title('y_2 d1');

subplot(6,2,2); plot(theta_21_a5); title('\theta_{21} a5');
subplot(6,2,4); plot(theta_21_d5); title('\theta_{21} d5');
subplot(6,2,6); plot(theta_21_d4); title('\theta_{21} d4');
subplot(6,2,8); plot(theta_21_d3); title('\theta_{21} d3');
subplot(6,2,10); plot(theta_21_d2); title('\theta_{21} d2');
subplot(6,2,12); plot(theta_21_d1); title('\theta_{21} d1');

% konstruiranje Wienerovog filtra
h_w_a5 = theta_21_a5.^2 ./ ( theta_21_a5.^2 + sigma_est^2 );
h_w_d5 = theta_21_d5.^2 ./ ( theta_21_d5.^2 + sigma_est^2 );
h_w_d4 = theta_21_d4.^2 ./ ( theta_21_d4.^2 + sigma_est^2 );
h_w_d3 = theta_21_d3.^2 ./ ( theta_21_d3.^2 + sigma_est^2 );
h_w_d2 = theta_21_d2.^2 ./ ( theta_21_d2.^2 + sigma_est^2 );
h_w_d1 = theta_21_d1.^2 ./ ( theta_21_d1.^2 + sigma_est^2 );

% prikaz Wienerovog filtra u domeni transformacije W2
figure('Name','Prikaz empirijskog Wienerovog filtra')
subplot(6,1,1); plot(h_w_a5); title('h_w a5');
subplot(6,1,2); plot(h_w_d5); title('h_w d5');
subplot(6,1,3); plot(h_w_d4); title('h_w d4');
subplot(6,1,4); plot(h_w_d3); title('h_w d3');
subplot(6,1,5); plot(h_w_d2); title('h_w d2');
subplot(6,1,6); plot(h_w_d1); title('h_w d1');

% filtracija Wienerovim filtrom u domeni transformacije W2
theta_2_a5 = y_2_a5.*h_w_a5;
theta_2_d5 = y_2_d5.*h_w_d5;
theta_2_d4 = y_2_d4.*h_w_d4;
theta_2_d3 = y_2_d3.*h_w_d3;
theta_2_d2 = y_2_d2.*h_w_d2;
theta_2_d1 = y_2_d1.*h_w_d1;

% inverzna transformacija W2

theta_2_C = [theta_2_a5   theta_2_d5   theta_2_d4   theta_2_d3   theta_2_d2   theta_2_d1];
theta_2_L = y_2_L;

% inverzna transformacija
s_est = waverec(theta_2_C, theta_2_L, W_2);

% prikaz rekonstruiranog signala nakon obrade Wienerovim filtrom
figure('Name','Prikaz signala s_est')
subplot(4,1,1), plot(s);
title('s'); xlabel('uzorak'); ylabel('amplituda');
subplot(4,1,2), plot(x);
title('x = s + n'); xlabel('uzorak'); ylabel('amplituda');
subplot(4,1,3), plot(s_est);
title('s est'); xlabel('uzorak'); ylabel('amplituda');
subplot(4,1,4), plot(s_est-s);
title('s est - s'); xlabel('uzorak'); ylabel('amplituda');








