% potrebuji pekne spektrum 

s = readl16('/homes/kazi/cernocky/SOUNDS/test.l16',0,inf); 
x = s(1350:(1350+511)); plot(x);
x = x .* hamming(512)'; plot(x);
X = fft(x); 
X = X(1:256); 
f = (0:255)/512 * 8000; 
plot (f, 10*log(abs(X).^2)); print -depsc FIG/spec_512.eps

% u cepstra by to chtelo postup ... 
x = s(1350:(1350+511)); plot(x);
plot(x); axis tight; print -depsc FIG_CEPS/signal.eps
cps = real(ifft(log(abs(fft(x)) .^ 2)));  % real potreba kvuli num. chybam 




% funkce na ceps Matlabu - odpovida tomu rceps. 
cpsm = rceps(x); 
plot (1:512, cps, 1:512, cpsm) % matlab je 2x mensi ???
plot (1:512, cps, 1:512, 2* cpsm) % matlab je 2x mensi ???
% ... ok, uplne stejny 
% pro print zoom ... 
plot(cps); print -depsc FIG_CEPS/ceps_src.eps

%%% trochu rozvest %%%
% signal uz mame. 
aux = abs(fft(x)) .^ 2
f = (0:255)/512 * 8000; plot (f, aux(1:256)); print -depsc FIG_CEPS/abspow2.eps
aux = log(abs(fft(x)) .^ 2); plot (f, aux(1:256)); print -depsc FIG_CEPS/log.eps
aux = real(ifft(log(abs(fft(x)) .^ 2))); plot (aux); axis tight;  print -depsc FIG_CEPS/ceps_all.eps
% vyrez mame v ceps. 

%%% hratky s nulovanim  - 31 vcetne c0, 30 vpravo !
cnofilter = cps; cnofilter(1:31) = 0; cnofilter((512-29):512) = 0; 
plot (cnofilter); axis tight; print -depsc FIG_CEPS/cnofilter.eps
% zpet do spektra real(ifft(log(abs(fft(x)) .^ 2)));
spnofilter = fft(cnofilter); plot (f, real(spnofilter(1:256))); 
  print -depsc FIG_CEPS/spnofilter.eps
% ted jeste signal ! budeme potrebovat faze z puvodniho spektra ... 
phase = angle (fft(x)); 
mag = sqrt(exp(real(spnofilter)));  plot (mag); print -depsc FIG_CEPS/magnofilter.eps
signofilter = real(ifft(mag .* exp(j* phase))); plot(signofilter); axis tight
  print -depsc FIG_CEPS/signofilter.eps
   
%% a ted naopak - bez buzeni - davame nulovou fazi - s puvodni to nehrne... 
cnoexc = cps; cnoexc(32:(512-30)) = 0; 
plot (cnoexc); axis tight;  print -depsc FIG_CEPS/cnoexc.eps
% zpet do spektra real(ifft(log(abs(fft(x)) .^ 2)));
spnoexc = fft(cnoexc); plot (f, real(spnoexc(1:256))); 
  print -depsc FIG_CEPS/spnoexc.eps
% ted jeste signal ! budeme potrebovat faze z puvodniho spektra ... 
phase = angle (fft(x)); 
mag = sqrt(exp(real(spnoexc)));  plot (mag); print -depsc FIG_CEPS/magnoexc.eps
signoexc = real(ifft(mag)); % .* exp(j* phase))); 
plot(signoexc); axis tight  % musim zoomovat ... 
  print -depsc FIG_CEPS/signoexc.eps