% zacatek uplne stejny jako minule ... 

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

% casovy odhad korelacni funkce - bereme 100. realizaci 
x = S(100,:); N = 320; 
R = xcorr(x,'biased'); k = -N+1:N-1
plot(k,R)
% ... ne, chce to delsi usek ... 
x = s(5000:5000+2047); N = 2048; 
R = xcorr(x,'biased'); k = -N+1:N-1
plot(k,R)
% ..., ne jedeme z jedne realizace ... 
x = S(100,:); N = 320; 
R = xcorr(x,'biased'); k = -N+1:N-1
plot(k,R); axis tight; grid
print -depsc FIG/biased.eps

R = xcorr(x,'unbiased'); k = -N+1:N-1
plot(k,R); axis tight; grid
print -depsc FIG/unbiased.eps

%%%%%%%%%%% psd jako dtft R(k) (viz aux.m z diskretni) 
Fs = 16000; 
f = -3*Fs:100:3*Fs; 
fn = f / Fs; 
Om = 2*pi*f;
om = Om / Fs; 
% dtft
R = xcorr(x,'biased'); k = -N+1:N-1
X = zeros(size(om)); 
for ii=1:length(om)
  ejon = exp(j*om(ii).*k); 
  X(ii) = sum(R .* ejon); 
end
subplot(211); plot (om,abs(X)); grid; axis tight; title('normalized omega');
subplot(212); plot (f,abs(X)); grid; axis tight;  title('f'); 
print -depsc FIG/psddtft.eps
% zoom (dost blbe, ale pujde to ...)
indi = find(om > -pi & om < pi); 
subplot(211); plot (om(indi),abs(X(indi))); grid; axis tight; title('normalized omega');
subplot(212); plot (f(indi),abs(X(indi))); grid; axis tight;  title('f'); 
% ... no to teda ne ... 

Fs = 16000; 
f = -Fs/2:Fs/512:Fs/2; 
fn = f / Fs; 
Om = 2*pi*f;
om = Om / Fs; 
% dtft
R = xcorr(x,'biased'); k = -N+1:N-1
X = zeros(size(om)); 
for ii=1:length(om)
  ejon = exp(j*om(ii).*k); 
  X(ii) = sum(R .* ejon); 
end
subplot(211); plot (om,abs(X)); grid; axis tight; title('normalized omega');
subplot(212); plot (f,abs(X)); grid; axis tight;  title('f'); 
print -depsc FIG/psddtftzoom.eps

%%% odhad pomoci DFT
x = S(100,:); N = 320; 
Gdft = 1/N * abs(fft(x)).^2; 
om = (0:159)/160 * pi; Gdft = Gdft(1:160); 
subplot(211); plot (om,Gdft); grid; title('one realization')

%%% averaging 
[Omega,Nfr] = size(S);    % je to blby, ale Nfr je delka frejmu ... 
bigg = zeros(Omega,160); 
for ii=1:Omega,
  x = S(ii,:);
  Gdft = 1/N * abs(fft(x)).^2; Gdft = Gdft(1:160); 
  bigg(ii,:) = Gdft; 
end
Gave = mean(bigg); 
subplot(212); plot (om,Gave); grid; title('averaging - Welch')
print -depsc FIG/psddft.eps 

%%%%%%%%%%%%%%%%% pruchod systemem %%%%%%%%%%%%%%%%%%%%%
x = S(100,:); N = 320; 
Gdft = 1/N * abs(fft(x)).^2; 
om = (0:159)/160 * pi; Gdft = Gdft(1:160); 
subplot(211); plot (x); grid; axis tight
subplot(212); plot (om,Gdft); grid; axis tight
print -depsc FIG/systin.eps 

% filtr - preem 
b = [1 -0.9]; a = 1; 
[H, om] = freqz(b,a,256,2*pi); 
subplot(211); plot(om,abs(H)); grid; axis tight; title('|H|'); 
subplot(212); plot(om,abs(H).^2); grid; axis tight; title('|H|^2'); 
print -depsc FIG/systh.eps 

% ted filtrovayt a pocitat psd vystupu. 
y = filter(b,a,x); N = 320; 
Gdft = 1/N * abs(fft(y)).^2; 
om = (0:159)/160 * pi; Gdft = Gdft(1:160); 
subplot(211); plot (y); grid; axis tight
subplot(212); plot (om,Gdft); grid; axis tight
print -depsc FIG/systout.eps 

%%%%%%%%%%%% gauss
x=-5:0.1:5;
px=zeros(size(x));
i=1;
for xx=x,
  px(i)=normal(xx, 0,1);
  i=i+1;
end
plot(x,px); title ('\mu=0, \sigma=1'); xlabel ('x'), ylabel ('p(x)'); grid
print -depsc FIG/gauss1.eps
Fx = cumsum(px) * 0.1; 
plot(x,Fx); title ('\mu=0, \sigma=1'); xlabel ('x'), ylabel ('F(x)'); grid
print -depsc FIG/Fgauss1.eps


x = randn(1,10000); 
plot(x(1:200)); 
print -depsc FIG/gwhite.eps

%% korelaky 
y = x(1:200); R = xcorr(y,'biased'); 
k = -199:199; plot(k,R); axis tight; grid; 
print -depsc FIG/Rwhite_bad.eps

%%% psd 
N = 200; 
Gdft = 1/N * abs(fft(y)).^2; 
om = (0:99)/100 * pi; Gdft = Gdft(1:100); 
subplot(211); 
plot (om,Gdft); grid; axis tight
%print -depsc FIG/Gwhite_1.eps 

%%% welch 
S = randn(1000,200); 
[Omega,Nfr] = size(S);    % je to blby, ale Nfr je delka frejmu ... 
bigg = zeros(Omega,100); 
for ii=1:Omega,
  x = S(ii,:);
  Gdft = 1/N * abs(fft(x)).^2; Gdft = Gdft(1:100); 
  bigg(ii,:) = Gdft; 
end
Gave = mean(bigg); 
subplot(212); plot (om,Gave); grid; title('averaging - Welch'); grid; axis tight
print -depsc FIG/Gwhite.eps 

%%%%%%%%%% quant %%%%%%%%%%%
x = 4 * cos(0.1 * (0:190)); 
subplot(311); plot(x); axis tight; grid; 
% q na 8 hodnot:
q = linspace(min(x),max(x),8)'; 
ind = q_code(x,q);
xq = q(ind);
subplot(312); plot(xq); axis tight; grid; 
% chyba
e = x - xq'; 
subplot(313); plot(e); axis tight; grid; 
print -depsc FIG/q.eps 

% pocitani snr (mel by byt 19.76 dB):
snr = 10*log10 (sum(x.^2) / sum(e.^2))

% ... 19.36 - ok. 
