% prvni cast neni matlab, ale nahravani signalu - 
% budu poustet vodu, potrebuju 1000 20ms useku (16 KHz) 
% to je 50 za sekundu, takze 20 sec. 

rec:
  -c, --channels=CHANNELS      specifies the number of sound channels in FILE
  -d, --device=DEVICE          use DEVICE for input/output
  -f, --format=FORMAT          specifies bit format of sample
                               FORMAT is either s, u, U, A, a, or g
  -r, --rate=RATE              sample rate in hertz of FILE
  -s, --size=SIZE              interpret size of sample
                               SIZE is either b, w, l, f, d, or D
  -t, --type=TYPE              specifies file format of FILE
  -v, --volume=VOLUME          change amplitude
  -x, --xinu                   reverse bit order of sample
                               (only works with 16-bit and 32-bit integer data)
      --file                   next argument is FILE
  -h, --help                   display this help and exit
      --version                output version information and exit
chci signed shorts ...
rec -f s -s w -r 16000 voda.wav 
mixer: mic 90% Igain 70% ... 
... hm, na konec jsem to nahral ve WIN a NA BATERKY (jinak je tam desny 
bordel). Mikrofon PRIMO NA VODOVODNI TRUBCE, jinak je to moc slabe. 

s = wavread('voda.wav'); 
% ramce 20 ms, jsme na 16k, takze 320. 
S = trame(s,320,0)'; % POZOR !!! CHCI JEDNOTLIVA MERENI V RADCICH ! 
[N,Nfr] = size(S);    % je to blby, ale Nfr je delka frejmu ... 
for ii=1:Nfr
  ii
  plot (S(ii,:)); axis([1 320 -1 1]); pause
end
% ... ok, vypadaji podobne ... 
% a few pictures for slides: 
t = (0:319) * 1/16000; n = 0:319; 
subplot(411); plot(t,S(1,:)); axis tight
subplot(412); plot(t,S(200,:)); axis tight
subplot(413); plot(t,S(500,:)); axis tight
subplot(414); plot(t,S(1000,:)); axis tight
print -depsc FIG/4real_spoj.eps
subplot(411); stem(n,S(1,:)); axis tight
subplot(412); stem(n,S(200,:)); axis tight
subplot(413); stem(n,S(500,:)); axis tight
subplot(414); stem(n,S(1000,:)); axis tight
print -depsc FIG/4real_disk.eps

%%%%%%% distrib fce a funkce hust. rozdeleni  pravd. 
% jak vypadaji histaky ? 
H = zeros(50,Nfr); 
[h,x] = hist(S(:,1),50); bar(x,h); 
for ii=1:Nfr
%  ii
  h = hist(S(:,ii),x); bar (x,h); %pause
  H(:,ii) = h';
end
% ... to je slusny, docela podobne ! 

%%%% par odhadu F(x) - 1 50 100 150
p = H(:,1)/N; F1 = cumsum(p);
p = H(:,50)/N; F50 = cumsum(p);
p = H(:,100)/N; F100 = cumsum(p);
p = H(:,150)/N; F150 = cumsum(p); 
plot (x,F1, x,F50, x,F100, x,F150); xlabel('x'); axis tight
legend('F(x,1)','F(x,50)','F(x,100)','F(x,150)'); grid
print -depsc FIG/F_disk.eps

[1/16000 50/16000 100/16000 150/16000]
% to stejne se spoj ... 
legend('F(x,0.1ms)','F(x,3.1ms)','F(x,6.3ms)','F(x,9.4ms)'); grid
print -depsc FIG/F_spoj.eps

%%% ted odhady p(x,t) ... ach jo, to zobrazovani ... 
Delta = x(2) - x(1)
L = 50; myx = [x-Delta/2; x+Delta/2]; myx = myx(:); 
p = H(:,1)/N/Delta; pzobr = [p'; p']; pzobr1 = pzobr(:); 
p = H(:,50)/N/Delta; pzobr = [p'; p']; pzobr50 = pzobr(:); 
p = H(:,100)/N/Delta; pzobr = [p'; p']; pzobr100 = pzobr(:); 
p = H(:,150)/N/Delta; pzobr = [p'; p']; pzobr150 = pzobr(:); 
subplot(221); plot(myx,pzobr1); xlabel('x'); ylabel('p(x,1)');axis tight
subplot(222); plot(myx,pzobr50); xlabel('x'); ylabel('p(x,50)');axis tight
subplot(223); plot(myx,pzobr100); xlabel('x'); ylabel('p(x,100)');axis tight
subplot(224); plot(myx,pzobr150); xlabel('x'); ylabel('p(x,150)');axis tight
print -depsc FIG/p_disk.eps

% jen prejmenovat ... 
subplot(221); plot(myx,pzobr1); xlabel('x'); ylabel('p(x,0.1ms)');axis tight
subplot(222); plot(myx,pzobr50); xlabel('x'); ylabel('p(x,3.1ms)');axis tight
subplot(223); plot(myx,pzobr100); xlabel('x'); ylabel('p(x,6.3ms)');axis tight
subplot(224); plot(myx,pzobr150); xlabel('x'); ylabel('p(x,9.4ms)');axis tight
print -depsc FIG/p_spoj.eps

% souborovy odhad strednich hodnot. 
mens = mean(S); plot (mens)
t = (0:319) * 1/16000; n = 0:319; 
subplot(211); plot(t,mens); grid; xlabel('t'); axis tight;
subplot(212); stem(n,mens); grid; xlabel('n'); axis tight;
print -depsc FIG/means.eps

% ... hm hm, zhruba stac ... 

% smer odchylek ... 
stds = std(S); 
subplot(211); plot(t,stds); grid; xlabel('t'); axis tight; v = axis;
subplot(212); stem(n,stds); grid; xlabel('n'); axis([0 319 v(3) v(4)]); 
print -depsc FIG/stds.eps
% ... hm hm, snad taky ... 


% a ted p(x1,x2,t1,t2) ... to jsou vlastne 2D histogramy, to bude maso !
% zvolime celkovy pocet bodu 10 x 10 = 100, aby se jeste vubec dalo neco 
% pocitat, interval -0.4 do 0.4 
% asi udelam funkci hist2 ... 

% napred zkusime pro cas 1 - 1
x = linspace(-0.4,0.4,30); 
h2 = hist2(S(:,1),S(:,1),x,x); 
imagesc (x,x,h2); axis xy; colorbar
% ... pekny, jen diag. 

% ted pro rozdil 0 az 20 vzorku, ale take pro cas 100 + 0 az 20 !
% pridam take odhad autokor. koeficientu podle teor vztahu 
% 5-8 ze Sebesty (numericka integrace). 
for k = 0:20,
  n1 = 1; n2=n1+k; n3=100; n4=n3+k; 
  [h1,p1,r1] = hist2(S(:,n1),S(:,n2),x,x); 
  [h2,p2,r2] = hist2(S(:,n3),S(:,n4),x,x); 
  %subplot(121); 
  imagesc (x,x,p1); axis xy; colorbar; xlabel('x2'); ylabel('x1');
  name = ['FIG/px1x2_' num2str(k) '.eps']; 
  print (1,'-depsc',name); 
%  subplot(122); imagesc (x,x,p2); axis xy; colorbar
  [k n1 n2 n3 n4 r1 r2]
  pause
end
% cool !!!! 

% compute it for k=0..40 and store somewhere. 
K = 40; 
R1 = zeros(1,K+1); R2 = zeros(1,K+1); 
for k = 0:K,
  n1 = 1; n2=n1+k; n3=100; n4=n3+k; 
  [h1,p1,r1] = hist2(S(:,n1),S(:,n2),x,x); 
  [h2,p2,r2] = hist2(S(:,n3),S(:,n4),x,x); 
  R1(k+1) = r1;   R2(k+1) = r2; k 
end
subplot(211); plot(0:K,R1); xlabel('k'); ylabel('R(0,k)'); grid
subplot(212); plot(0:K,R2); xlabel('k'); ylabel('R(100,100+k)'); grid 
print -depsc FIG/R0R100.eps

% compare with estimates from time domain (use whole signal)

Rtim = xcorr(s,'biased');