close all;
clear;
clc;

fs=8000;
M=50;
fc = fs/4;
Wc = 0.5;

h_lp = fir1(M, Wc, 'low');
H_len_half = ceil(length(h_lp)/2);

if mod(H_len_half, 2) == 0
  mask = [-1;1];
else
  mask = [1;-1];
end

mask = repmat(mask, H_len_half);
mask = mask(1:length(h_lp));

h_hp = h_lp .* mask;

a=1;
b=h_hp;

figure();
freqz(b,a);
figure();
impz(b,a);

b_hp = h_hp;
b_lp = h_lp;

sig=loadbin("s0001.bin");

% HP filt
y_hp=filter(b,a,sig);
figure();
subplot(211);
spectrogram(y_hp,256,[],[],'yaxis');
colormap jet;
title("HP");
% DP filt
y_lp=filter(h_lp,a,sig);
subplot(212);
spectrogram(y_lp,256,[],[],'yaxis');
colormap jet;
title("LP");

y_rekonst_2=y_lp+y_hp; % secteni filtrovanych signalu pro dalsi ucely

y_dec_hp=y_hp(1:2:length(y_hp));
y_dec_lp=y_lp(1:2:length(y_lp));

% HP decim. plot
figure();
subplot(211);
spectrogram(y_dec_hp,128,[],[],'yaxis');
colormap jet;
title("HP decimated");
% DP decim. plot
subplot(212);
spectrogram(y_dec_lp,128,[],[],'yaxis');
colormap jet;
title("LP decimated");

% naplneni nulami
y_dec_hp_0s = zeros(1, 2*length(y_dec_hp));
y_dec_lp_0s = zeros(1, 2*length(y_dec_lp));
y_dec_hp_0s(1:2:end) = y_dec_hp;
y_dec_lp_0s(1:2:end) = y_dec_lp;

% HP decim. plot
figure();
subplot(211);
spectrogram(y_dec_hp_0s,256,[],[],'yaxis');
colormap jet;
title("interpol. HP");
% DP decim. plot
subplot(212);
spectrogram(y_dec_lp_0s,256,[],[],'yaxis');
colormap jet;
title("LP");

% filtrace
y_dec_hp_filt=filter(b_hp,a,y_dec_hp_0s);
y_dec_lp_filt=filter(b_lp,a,y_dec_lp_0s);
y_rekonst=2*(y_dec_hp_filt+y_dec_lp_filt);

figure();
subplot(211);
spectrogram(y_rekonst,256,[],[],'yaxis');
colormap jet;
title("Rekonstruovany");
subplot(212);
spectrogram(sig,256,[],[],'yaxis');
colormap jet;
title("puvodni");

figure();
plot(sig);
hold on;
plot(y_rekonst);
hold off;

% vypocet SNR
figure();
posuv=M;
sig=sig(:);
y_rekonst=y_rekonst(:);
sig_a=sig(1:end-posuv);
sig_b=y_rekonst(posuv+1:end-1);
error=sig_b-sig_a;
plot(error);
pow_a = mean(sig_a.^2);
pow_error = mean(error.^2);

snr=10*log10(pow_a/pow_error)

% vypocet signalu bez decimace a interpolace
posuv = M/2;
sig_a=sig(1:end-posuv);
sig_b=y_rekonst_2(posuv+1:end);
error=sig_b-sig_a;
plot(error);
pow_a = mean(sig_a.^2);
pow_error = mean(error.^2);
snr2=10*log10(pow_a/pow_error)