clear
close all

%fs = 44.1e3;
%N = 120*fs;
%in = conv( randn(N,1),exp(-5*(1e-2)*(1:100)'), 'same' );

[in, fs] = audioread('Sounds/master.mp3');
% Master of puppets
in =  sum( in(235*fs:248*fs,:) , 2 );
in = (in+0.1*randn(length(in),1)).*tukeywin(length(in),0.175);
N = length(in);
%%
Na = 256*8;
in = in(1:floor(length(in)/Na)*Na);

r = 0;
Nstep = floor ( N / Na );

for n = 1 : Nstep
    win = hann(Na).*in( (n-1) * Na + 1 : n * Na );
    r = r + xcorr( win, win, 'coeff'  )/Nstep;
end

r0 = r(floor(length(r)/2) + 2 : floor(length(r)/2) + Na  );
R0 = r(floor(length(r)/2) + 1 : floor(length(r)/2) + Na - 1  );
R = zeros( Na-1, Na-1 );
for n = 1 : Na-1
    R( n, : ) = r( floor(length(r)/2) + 2 - n  : floor(length(r)/2) + Na - n ) ;
end
%
a = R\r0;
%encoded_out = zeros(N,1);
% for n = Na : N - 1
%   encoded_out(n) = in(n) - flipud(a).'*in( n - Na + 1 : n  - 1 );
% end
encoded_out = filter([1;-a],1,in);
%%
h = [1;-a];

figure('units','normalized','outerposition',[0 0 1 1])
subplot(1,2,1);
plot((0:Na-1)/fs,h,'LineWidth',1.5);
xlabel('t -> [s]');
ylabel('h -> [V]');
xlim([-2/fs,(Na-2)/fs*0.2])
title('Ideális predikciós szûrõ impulzusválasza')
grid on

subplot(1,2,2);
semilogx((0:Na-1)/Na*fs,20*log10(abs(fft(h))),'LineWidth',1.5);
xlabel('f -> [Hz]');
ylabel('|H| -> [dB]');
xlim([0,fs/2])
title('Predikciós szûrõ átviteli függvénye')
grid on
%%

in_mx = reshape(in,Na,length(in)/Na);
in_sp = sum(abs(fft(in_mx,[],1)),2)/Nstep;

out_mx = reshape(encoded_out,Na,length(encoded_out)/Na);
out_sp = sum(abs(fft(out_mx,[],1)),2)/Nstep;


r_in = 0;
r_out = 0;
for n = 1 : Nstep
    win1 = hann(Na).*in( (n-1) * Na + 1 : n * Na );
    win2 = hann(Na).*encoded_out( (n-1) * Na + 1 : n * Na );
    r_in = r_in + xcorr( win1, win1, 'coeff'  )/Nstep;
    r_out = r_out + xcorr( win2, win2, 'coeff'  )/Nstep;
end
%%
figure('units','normalized','outerposition',[0 0 1 1])
subplot(2,2,1);
plot((1:length(in))/fs,in,'LineWidth',1.5);
hold on
plot((1:length(in))/fs,encoded_out,'LineWidth',1.5);
xlabel('t -> [sec]');
ylabel('A -> [V]');
grid on
legend('Original signal','Residual')
title('Time histories')

subplot(2,2,2);
semilogx((0:Na-1)/Na*fs,20*log10(abs(fft(in_sp))),'LineWidth',1.5);
hold on
semilogx((0:Na-1)/Na*fs,20*log10(abs(fft(out_sp))),'LineWidth',1.5);
xlabel('f -> [Hz]');
ylabel('|A| -> [dB]');
xlim([0,fs/2])
grid on
legend('Original signal','Resiudal');
title('Average spectra')

dt = ((1:(2*Na-1))-round((2*Na-1)/2))/(2*Na-1)/fs;
subplot(2,2,3);
plot(dt,r_in,'LineWidth',1.5);
hold on
plot(dt,r_out,'LineWidth',1.5);
xlabel('dt -> [sec]');
ylabel('A -> [V]');
xlim([dt(1),dt(end)])
grid on
legend('Original signal','Residual')
title('Autocorrelation of signals')

subplot(2,2,4);
semilogx((1:2*Na-1)/(2*Na-1)*fs,  20*log10(abs(fft(r_in))),'LineWidth',1.5);
hold on
semilogx((1:2*Na-1)/(2*Na-1)*fs,20*log10(abs(fft(r_out))),'LineWidth',1.5);
xlabel('f -> [Hz]');
ylabel('|A| -> [dB]');
xlim([100,fs/2])
grid on
legend('Original signal spectrum','Resiudal spectrum');
title('Spectral power densities')


%% Quantization

q = 3;
Q = 10;
encoded_out_2 = encoded_out*Q;
encoded_out_2(abs(encoded_out_2)>1) = 0;
encoded_out_q = round( encoded_out_2*(2^q-1)) / (2^q-1)/Q;
in = in/max(in);
in_q = round( in*(2^q-1) ) / (2^q-1);


E = encoded_out - encoded_out_q;

plot(filter(1,[1;-a],E))
%
Hd = dfilt.df1(1,[1;-a]);

%%
% Decoder
%decoded_out = zeros(N,1);
%for n = Na : N - 1
%  decoded_out(n) = encoded_out_q(n) + flipud(a).'*decoded_out( n - Na + 1 : n - 1 );
%end
decoded_out= filter(1,[1;-a],encoded_out_q);

figure
plot(  running_rms ( decoded_out(1:length(in) - in) , fs, 0.1 ));