%%%%%%%% dft obd impulsu %%%%%%%%%%%
N = 16; 
n = 0:N-1; 
x = zeros(size(n)); x(1:4) = 1; x(13:16)=1;  
subplot(311); stem(n,x); 
X = fft(x); k = n; 
subplot(312); stem(k,abs(X)); grid
subplot(313); stem(k,angle(X)); grid
print -depsc FIG/dftobd.eps

% ruzne frekv osy pro Fs = 44100 Hz: 
Fs = 44100; 
fn = k/N; 
f = fn * Fs; 
omn = 2*pi*fn; 
om = omn * Fs; 
subplot(221); stem(fn,abs(X)); grid; axis([0 1 0 8]); title('norm. f'); 
subplot(222); stem(f,abs(X)); grid;  axis([0 Fs 0 8]); title('f');
subplot(223); stem(omn,abs(X)); grid; axis([0 2*pi 0 8]); title('norm. omega');
subplot(224); stem(om,abs(X)); grid;  axis([0 2*pi*Fs 0 8]); title('omega');
print -depsc FIG/fosy.eps

%%%%%%%%%%%% 1 per harmon 
N = 16; 
n = 0:N-1; 
x = cos(2*pi/N * n); 
subplot(311); stem(n,x); 
X = fft(x); k = n; 
subplot(312); stem(k,abs(X)); grid
%ang = angle(X) .* sign(X); % zeroing unwanted stuff ... does not work
ang = zeros(size(k));
subplot(313); stem(k,ang); grid
print -depsc FIG/dftharm.eps

% ruzne frekv osy pro Fs = 44100 Hz: 
Fs = 44100; 
fn = k/N; 
f = fn * Fs; 
omn = 2*pi*fn; 
om = omn * Fs; 
subplot(221); stem(fn,abs(X)); grid; axis([0 1 0 8]); title('norm. f'); 
subplot(222); stem(f,abs(X)); grid;  axis([0 Fs 0 8]); title('f');
subplot(223); stem(omn,abs(X)); grid; axis([0 2*pi 0 8]); title('norm. omega');
subplot(224); stem(om,abs(X)); grid;  axis([0 2*pi*Fs 0 8]); title('omega');
print -depsc FIG/fosyharm.eps

%%%%%%%%%%%%%%%%%%%% cykl posun N = 16; 
n = 0:N-1; 
x = cos(2*pi/N * n); 
subplot(321); stem(n,x); 
X = fft(x); k = n; 
subplot(323); stem(k,abs(X)); grid
%ang = angle(X) .* sign(X); % zeroing unwanted stuff ... does not work
ang = zeros(size(k));
subplot(325); stem(k,ang); grid; axis([0 N-1 -pi pi])
y = cos(2*pi/N * (n-2)); 
subplot(322); stem(n,y); 
Y = fft(y); k = n; 
subplot(324); stem(k,abs(Y)); grid
%% needs some processing !
ang = zeros(size(k)); ang(2) = angle(Y(2)); ang(16) = angle(Y(16)); 
subplot(326); stem(k,ang); grid; axis([0 N-1 -pi pi])
print -depsc FIG/cyklpos.eps
% kontrola pro 1: 
2*pi*1*2 / 16 
% ... 0.7854 ... ok. 

%%%%%%%%%%%%%%%%%%%% kruh konv. 
n = 0:3;
x1 = [2 2 0 0]; 
x2 = [1 -1 0 0]; 
X1 = fft(x1); X2 = fft(x2); 
subplot(331); stem (n,x1); grid
subplot(334); stem (n,abs(X1)); grid
subplot(337); stem (n,angle(X1)); grid
subplot(332); stem (n,x2); grid
subplot(335); stem (n,abs(X2)); grid
subplot(338); stem (n,angle(X2));  grid
X = X1 .* X2; x = ifft(X); 
subplot(333); stem (n,real(x)); grid
subplot(336); stem (n,abs(X)); grid
% again hack for the phase ... subplot(339); stem (n,angle(X));  grid
subplot(339); stem (n,zeros(1,4));  grid
print -depsc FIG/kruhkonv.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% analog
N = 16; n = 0:15;  T = 0.001; nT = n*T;  t = 0:0.0001:16*0.001; 
omn = 2*pi/N; om = omn / T; 
x = 10*cos(om * t + pi/4); xn = 10*cos(omn * n + pi/4); 
subplot(311); plot(t,x); hold on; stem(nT,xn,'r'); hold off; grid
Xk = fft(xn); 
subplot(312); stem(n,abs(Xk)); grid
ang = zeros(1,N); ang(2) = angle(Xk(2)); ang(16) = angle(Xk(16)); 
subplot(313); stem(n,ang); grid
print -depsc FIG/harm+cont.eps

% analog - bad spatna perioda 
N = 16; n = 0:15;  T = 0.001; nT = n*T;  t = 0:0.0001:16*0.001; 
omn = 1.2*2*pi/N; om = omn / T; 
x = 10*cos(om * t); xn = 10*cos(omn * n); 
subplot(311); plot(t,x); hold on; stem(nT,xn,'r'); hold off; grid
Xk = fft(xn); 
subplot(312); stem(n,abs(Xk)); grid
subplot(313); stem(n,angle(Xk)); grid
print -depsc FIG/badharm+cont.eps

% analog obdelnik 
N = 64; n = 0:N-1;  T = 0.001; nT = n*T;  t = 0:0.0001:N*0.001; 
theta = 0.032; T1 = N*T;  
x = zeros(1,length(t)); 
x(find(t >= T1/2 - theta/2 & t <= T1/2 + theta/2)) =1;
xn = zeros(1,N); 
xn(find(n >= N/2 - 16 & n <= N/2 + 16)) =1;
omn = 2*pi/N; om = omn / T; 
subplot(311); plot(t,x); hold on; stem(nT,xn,'r'); hold off; grid; axis tight
Xk = fft(xn); 
subplot(312); stem(n,abs(Xk)); grid; axis tight
subplot(313); stem(n,angle(Xk)); grid; axis tight
print -depsc FIG/obd+cont.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% analog obdelnik 
N = 64; n = 0:N-1;  T = 0.001; nT = n*T;  t = 0:0.0001:N*0.001; 
theta = 0.032; T1 = N*T;  
x = zeros(1,length(t)); 
x(find(t >= T1/2 - theta/2 & t <= T1/2 + theta/2)) =1;
% puvodni od minus 32ms
subplot(311); plot(t-0.032,x); grid; axis tight
xn = zeros(1,N); 
xn(find(n >= N/2 - 16 & n <= N/2 + 16)) =1;
omn = 2*pi/N; om = omn / T; 
subplot(312); plot(t,x); grid; axis tight
subplot(313); plot(t,x); hold on; stem(nT,xn,'r'); hold off; grid; axis tight
print -depsc FIG/obd+cont+posun.eps

% spek 
npul = 1:N/2; 
Xk = fft(xn); 
subplot(211); plot(n(npul),abs(Xk(npul))); grid; axis tight
subplot(212); plot(n(npul),angle(Xk(npul))); grid; axis tight
print -depsc FIG/ftobd.eps

% zero pad 
xpad = [xn zeros(1,512-64)];  subplot(311);  stem(xpad); 
Xk = fft(xpad); 
k = 0:255; ind = 1:256; 
subplot(312); plot(k,abs(Xk(ind))); grid; axis tight
subplot(313); plot(k,angle(Xk(ind))); grid; axis tight
print -depsc FIG/ftobd_pad.eps

% slusna fr. osa a oprava faze, velikost
om = 2*pi*k/512/T; 
Xopr = T * Xk(ind) .* exp(j * 0.032 * om); 
subplot(211); plot(om,abs(Xopr)); grid; axis tight
subplot(212); plot(om,angle(Xopr)); grid; axis tight
print -depsc FIG/frosa+fazeopr.eps




%%%