% pitch - rpiprvaa 2004 - prevzato ze stare, neco dodelane,
% vetsinou ne print ! - hm, asi to udelam trochu novy ... na konci 
%  souboryu nechavam, co bylo predtim ! 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dlouhodoby prediktor . 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s = readl16('/homes/kazi/cernocky/SOUNDS/test.l16', 0, inf); 
s = s(50:length(s)); 
s = s- mean(s); axis tight; 
plot (s); print -depsc FIG/cely_sig.eps

% zobrazeni ACF (delsi, tak aby tam bylo nekolik lagu)
odkud=1100; N=320;
x = s(odkud:odkud+N-1); plot (x); axis tight; 
plot (x); print -depsc FIG/frame.eps

% ACF 
cc = xcorr (x); 
cc = cc (N:length(cc)); m = 0:N-1; 
plot (m,cc); xlabel ('m'); ylabel ('R(m)');  axis tight; 
print -deps FIG/acf.eps

% ilustrace posunuti 
pos = [zeros(1,87) x]; orig = [x zeros(1,87)]; 
subplot(211); plot (1:(320+87),orig); axis tight; 
subplot(212); plot (1:(320+87),pos);  axis tight; 
print -deps FIG/posun.eps

%%% znazroneni lagu - uprava acf.eps v xfigu => acfmax.eps

% amdf:
amdf = zeros(1,319);
for nn = 1:319,
  posod = nn; posdo = 320; len = posdo - posod +1
  ood = 1; odo = 1 + len -1; 
  amdf(nn) = sum(abs( x(posod:posdo) - x(ood:odo))); 
end
plot (amdf); axis tight; print -deps FIG/amdf.eps

%%% NCCF - ukazka posunu - musim domalovat meze !
odkud=1100; N=320;
orig  = s(odkud:odkud+N-1+87); 
pos   = s(odkud-87:odkud+N-1);
subplot(211); plot (1:(320+87),orig); axis tight; 
subplot(212); plot (1:(320+87),pos);  axis tight; 
print -deps FIG/nccf_src.eps

% => nccf_pos.eps

% malovat energie nebudu ... anebo jo ? 
odkud=9800; N=320;
orig  = s(odkud:odkud+N-1+120); 
pos   = s(odkud-120:odkud+N-1);
subplot(211); plot (1:(320+120),orig); axis ([1,320+120,-1500,3000]); 
subplot(212); plot (1:(320+120),pos);  axis ([1,320+120,-1500,3000]); 
print -deps FIG/bad_src.eps

% => bad_case.eps

% zase pro puvodni pomoci NCCF  - good case
odkud=1100; N=320;
ccf = zeros(1,320); 
nccf = zeros(1,320); 
for k=0:319
  orig=s(odkud:odkud+N-1);
  pos=s(odkud-k:odkud+N-1-k);
  ccf(k+1) = sum(orig .* pos); 
  e1 = sum(orig.^2);   e2 = sum(pos.^2); 
  nccf(k+1) = ccf(k+1) / sqrt(e1 * e2);
end
subplot(211); plot (ccf); axis tight; 
subplot(212); plot (nccf); axis tight; 
print -depsc FIG/ccf_nccf_ok.eps

% zase pro puvodni pomoci NCCF  - good case
odkud=9800; N=320;
ccf = zeros(1,320); 
nccf = zeros(1,320); 
for k=0:319
  orig=s(odkud:odkud+N-1);
  pos=s(odkud-k:odkud+N-1-k);
  ccf(k+1) = sum(orig .* pos); 
  e1 = sum(orig.^2);   e2 = sum(pos.^2); 
  nccf(k+1) = ccf(k+1) / sqrt(e1 * e2);
end
subplot(211); plot (ccf); axis tight; 
subplot(212); plot (nccf); axis tight; 
print -depsc FIG/ccf_nccf_bad.eps

% to by snad stacilo ... 


%%%%%%%%%%%%%%%% +++ bylo +++ %%%%%%%%%%%% 

s = readl16('/homes/kazi/cernocky/SOUNDS/test.l16', 0, inf); 
s = s(50:length(s)); 
s = s- mean(s); 

odkud = 1800; N = 160;
x = s(odkud:odkud+N-1);
figure(1); plot (x); xlabel ('samples'); print -deps FIG/orig.eps

a = real(lpc (x, 10)); 
e = filter (a,1,x);
figure (2); plot(e); xlabel ('samples'); print -deps FIG/short.eps
ei = interp (e,6)

% korelace chyby ... s interpolaci 
%cc = xcorr (ei); 
%cc = cc (6*N:length(cc));
%figure (3); plot (cc); 

% bez interpolace 
cc = xcorr (e); 
cc = cc (N:length(cc));
figure (3); plot (cc); 

%%%%%%%%%% chce to pred. 2. radu !!!!!!!!!! %%%%%%%%%
%jedeme podle funkce K(m)
% NEbudeme upsamplovat ! 
%close all; 

% pocitame K(m) podle vztahu K(m) = r^2(m) + r^2(m-1) - 2 r(1) r (m-1) r (m);
r = cc / cc(1); 
rposun = [0 r(1:length(r)-1)]; 
K = r.^2 + rposun.^2 - 2 * r(2) * r .* rposun; 
figure (4); plot (K) 

% vymazeme prvnich 10, najdeme max a prevedeme jej na lag, 
% v m nechavame Matlabovsky index !
K(1:10) = zeros(1,10); 
[Km,m] = max (K); 
L=m-1

% udelej koeficienty filtru a vlastni filtr:
b1 = (r(2)*r(m) - r(m-1))/ (1 - r(2)^2)
b2 = (r(2)*r(m-1) - r(m))/ (1 - r(2)^2)
B = [1 zeros(1,L-2) b1 b2];

% filtruj residual: 
ee = filter (B,1,e); 
figure(5); plot (ee); xlabel ('samples'); print -deps FIG/long.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dalsi triky: 

% zobrazeni ACF (delsi, tak aby tam bylo nekolik lagu)
odkud=1100; N=320;
x = s(odkud:odkud+N-1); plot (x);
cc = xcorr (x); 
cc = cc (N:length(cc)); m = 0:N-1; 
figure (3); plot (m,cc); xlabel ('m'); ylabel ('R(m)'); 
print -deps FIG/Rm.eps

%%%%%%%%%%%%%%
% klipovani: na kratkem  framu: 
close all
%x = x (1:160); 
figure(1); plot (x);  xlabel ('samples'); ylabel ('original');
print -deps FIG/noclip.eps

cl = 0.6 * max (abs (x)); 
x1 = x;
x1 (find(x1 <= cl & x1 >= -cl)) = 0; 
kladny = find (x1 > 0 ); 
x1 (kladny) = x1 (kladny) - cl; 
zaporny = find (x1 > 0 ); 
x1 (zaporny) = x1 (zaporny) + cl; 
figure (2); plot (x1);  xlabel ('samples'); ylabel ('clip 1');
print -deps FIG/clip1.eps

x2 = sign (x1); 
figure (3); plot (x2);  xlabel ('samples'); ylabel ('clip 2');
v = axis; axis ([v(1) v(2) -1.1 1.1])
print -deps FIG/clip2.eps

%%%%%% korelace klipacu 
cc1 = xcorr (x1); 
cc1 = cc1 (N:length(cc1)); m = 0:N-1; 
figure (3); plot (m,cc1); xlabel ('m'); ylabel ('Rcl1(m)'); 
print -deps FIG/Rcl1.eps

cc2 = xcorr (x2); 
cc2 = cc2 (N:length(cc2)); m = 0:N-1; 
figure (3); plot (m,cc2); xlabel ('m'); ylabel ('Rcl2(m)'); 
% rcl2 je uplne stejny , takze ho netiskeneme ! print -deps FIG/Rcl2.eps

%%%%% korelace chyboveho 
a = lpc (x, 10); 
e = filter (a,1,x);
cce = xcorr (e); 
cce = cce (N:length(cce)); m = 0:N-1; 
figure (3); plot (m,cce); xlabel ('m'); ylabel ('Re(m)'); 
print -deps FIG/Re.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cps 
c = real (ifft (log (abs(fft(x)).^2))); % zasrany num chyby ...
c = c (1:N/2);
plot (c(20:160)) 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ted urceni pro cely signal, uf ... 
FR = trame (s, 240, 80); 
[N, Nfr] = size (FR)
Lall = []; rall = []; 
for i = 1:Nfr, 
  x = FR(:,i);
  [L,R] = pitchcor (x, 20,159, 0.3); 
  r = R / R(1); 
  Lall = [Lall L]; 
  rall = [rall r];
end
 plot (Lall); 

% troska zobrazeni : 
rall (rall < 0 ) = 0; % minusy nas netankuji
rall (1,:) = 0;    % toto nas nezajima 
imagesc (-rall); colormap('gray'); axis xy; xlabel('frame'); ylabel ('m');
title ('autocorrelation functions in all frames')

print -deps FIG/rall.eps 


% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
% ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

dal uz neprovadet !!!

%%%%% primo z korel koef -- nic moc !!!

% pred. koeficient a lag se asi pocitaji takhle:
R0 = cc(1); cc(1:10) = zeros(1,10); 
[Rm,L] = max (cc); 
L=L-1 
beta = - Rm/R0

% zkusime filtrovat: 
ee = filter ([1 zeros(1,L-1) beta],1,e);
figure(4); plot (ee);