%ukazka aliasingu - jo nebo nee. 

Fs = 8000; fmax = 3000; 
Os = 16000 * pi; omax = 6000 * pi; 
o = (-400:400)/100 * Os; % presne 100 bodu na jednu periodu. 
gj0 = find (o==0); 

zind = (-300:300) + gj0;  %indexy na zobr. 
zob = o(zind);   % omegy na zobr. 

%%%%%%%%%%%%%%% prvni priklad - no aliasing %%%%%%%%
%ind = (find(o==-omax):find(o==omax)); ... hm, tohle nejde (mat, presnost, tak rucne ...
X = zeros(size(o)); 
kam = (-37:37) + gj0; X(kam) = hamming(2*37+1); 
subplot(411); plot (zob, X(zind)); axis tight
% ted periodizovane - okolo -3,-2,-1,0,1,2,3
T = 1/Fs; 
Xs1 = zeros(size(o)); Xs2=Xs1; Xs3=Xs1;Xs4=Xs1;Xs5=Xs1;Xs6=Xs1; Xs7=Xs1;  
kam = (-37:37) + gj0 - 300; Xs1(kam) = 1/T * hamming(2*37+1);
kam = (-37:37) + gj0 - 200; Xs2(kam) = 1/T * hamming(2*37+1);
kam = (-37:37) + gj0 - 100; Xs3(kam) = 1/T * hamming(2*37+1);
kam = (-37:37) + gj0 - 0; Xs4(kam) = 1/T * hamming(2*37+1);
kam = (-37:37) + gj0 + 100; Xs5(kam) = 1/T * hamming(2*37+1);
kam = (-37:37) + gj0 + 200; Xs6(kam) = 1/T * hamming(2*37+1);
kam = (-37:37) + gj0 + 300; Xs7(kam) = 1/T * hamming(2*37+1);
Xs  = Xs1 + Xs2 + Xs3 + Xs4 + Xs5 + Xs6 + Xs7; 
subplot(412); plot(zob,Xs1(zind),'b',zob,Xs2(zind),'b',zob,Xs3(zind),'b',zob,Xs4(zind),'b',zob,Xs5(zind),'b',zob,Xs6(zind),'b',zob,Xs7(zind),'b',zob,Xs(zind),'r'); axis tight
% filtr
H = zeros(size(o));
H((-50:50) + gj0) = T; 
subplot(413); plot(zob,H(zind)); axis tight
%
Xr = H .* Xs; 
subplot(414); plot(zob,Xr(zind)); axis tight
print -depsc FIG/pr_noalias.eps

%%%%%% ted s aliasem  -- uz ne 3000, ale 7000 Hz ! 
X = zeros(size(o)); 
kam = (-87:87) + gj0; X(kam) = hamming(2*87+1); 
subplot(411); plot (zob, X(zind)); axis tight
% ted periodizovane - okolo -3,-2,-1,0,1,2,3
T = 1/Fs; 
Xs1 = zeros(size(o)); Xs2=Xs1; Xs3=Xs1;Xs4=Xs1;Xs5=Xs1;Xs6=Xs1; Xs7=Xs1;  
kam = (-87:87) + gj0 - 300; Xs1(kam) = 1/T * hamming(2*87+1);
kam = (-87:87) + gj0 - 200; Xs2(kam) = 1/T * hamming(2*87+1);
kam = (-87:87) + gj0 - 100; Xs3(kam) = 1/T * hamming(2*87+1);
kam = (-87:87) + gj0 - 0; Xs4(kam) = 1/T * hamming(2*87+1);
kam = (-87:87) + gj0 + 100; Xs5(kam) = 1/T * hamming(2*87+1);
kam = (-87:87) + gj0 + 200; Xs6(kam) = 1/T * hamming(2*87+1);
kam = (-87:87) + gj0 + 300; Xs7(kam) = 1/T * hamming(2*87+1);
Xs  = Xs1 + Xs2 + Xs3 + Xs4 + Xs5 + Xs6 + Xs7; 
subplot(412); plot(zob,Xs1(zind),'b',zob,Xs2(zind),'b',zob,Xs3(zind),'b',zob,Xs4(zind),'b',zob,Xs5(zind),'b',zob,Xs6(zind),'b',zob,Xs7(zind),'b',zob,Xs(zind),'r'); axis tight
% filtr
H = zeros(size(o));
H((-50:50) + gj0) = T; 
subplot(413); plot(zob,H(zind)); axis tight
%
Xr = H .* Xs; 
subplot(414); plot(zob,Xr(zind)); axis tight
print -depsc FIG/pr_alias.eps

%%%%%%%%%%%%% aa filtr %%%%%%%%%%%%%%%%%%%%%%
X = zeros(size(o)); 
kam = (-87:87) + gj0; X(kam) = hamming(2*87+1); 
subplot(511); plot (zob, X(zind)); axis tight
% aa
aux = hamming(2*87+1); aux = aux (88-49:88+49); 
Xaa = zeros(size(o)); 
kam = (-49:49) + gj0; Xaa(kam) = aux; 
subplot(512); plot (zob, Xaa(zind)); axis tight

% ted periodizovane - okolo -3,-2,-1,0,1,2,3
T = 1/Fs; 
Xs1 = zeros(size(o)); Xs2=Xs1; Xs3=Xs1;Xs4=Xs1;Xs5=Xs1;Xs6=Xs1; Xs7=Xs1;  
kam = (-49:49) + gj0 - 300; Xs1(kam) = 1/T .* aux;
kam = (-49:49) + gj0 - 200; Xs2(kam) = 1/T .* aux;
kam = (-49:49) + gj0 - 100; Xs3(kam) = 1/T .* aux;
kam = (-49:49) + gj0 - 0; Xs4(kam) = 1/T .* aux;
kam = (-49:49) + gj0 + 100; Xs5(kam) = 1/T .* aux;
kam = (-49:49) + gj0 + 200; Xs6(kam) = 1/T .* aux;
kam = (-49:49) + gj0 + 300; Xs7(kam) = 1/T .* aux;
Xs  = Xs1 + Xs2 + Xs3 + Xs4 + Xs5 + Xs6 + Xs7; 
subplot(513); plot(zob,Xs1(zind),'b',zob,Xs2(zind),'b',zob,Xs3(zind),'b',zob,Xs4(zind),'b',zob,Xs5(zind),'b',zob,Xs6(zind),'b',zob,Xs7(zind),'b',zob,Xs(zind),'r'); axis tight
% filtr
H = zeros(size(o));
H((-50:50) + gj0) = T; 
subplot(514); plot(zob,H(zind)); axis tight
%
Xr = H .* Xs; 
subplot(515); plot(zob,Xr(zind)); axis tight
print -depsc FIG/pr_aa.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reko - cas 
T = 1/8000; Fs = 8000; Os = 2*pi * Fs; 
t = -10*T:T/10:20*T;
zind = find (t>=0 & t<=10*T);  % pro zobrazeni jen vysek, aby to bylo pekne !
zt= t(zind); 
% signal 
x = sin(2*pi*600*t); 
subplot(411); plot (zt,x(zind)); axis tight; 
% vzorkovany
n = -10:20; nT = n*T;    zindn=11:21; % indexy na zobrazeni !  
xn = sin(2*pi*600*nT);   
subplot(412); stem(nT(zindn), xn(zindn)); axis tight;
% ted vygenerujeme mraky interp. signalu: 
HR = zeros(31,length(t)); 
for ii=n, 
  nii = ii+11; 
  hr = xn(nii)*sinc(Os/2 * (t-ii*T) / pi); HR(nii,:) = hr; % debilni definice sinc ! 
end
% priklad jednoho
subplot(413);
plot (zt,HR(14,zind),'b'); axis tight; hold on; stem(nT(14),xn(14)); hold off
% a ted vse, plus suma: 
xr= sum(HR); 
subplot(414);
plot (zt,HR(:,zind),'b',zt,xr(zind),'r'); axis tight; hold on; 
stem(nT(zindn),xn(zindn)); hold off
print -depsc FIG/reko.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%% ted pro obd. impuls od 2T do 5T %%%
T = 1/8000; Fs = 8000; Os = 2*pi * Fs; 
t = -10*T:T/10:20*T;
zind = find (t>=0 & t<=10*T);  % pro zobrazeni jen vysek, aby to bylo pekne !
zt= t(zind); 
% signal 
x = zeros(size(t)); x(find(t>=2*T & t<=5*T)) = 1; 
subplot(411); plot (zt,x(zind)); axis tight; 
% vzorkovany
n = -10:20; nT = n*T;    zindn=11:21; % indexy na zobrazeni !  
xn = zeros(size(n)); xn((2:5) + 11) = 1; 
subplot(412); stem(nT(zindn), xn(zindn));  axis tight;
% ted vygenerujeme mraky interp. signalu: 
HR = zeros(31,length(t)); 
for ii=n, 
  nii = ii+11; 
  hr = xn(nii)*sinc(Os/2 * (t-ii*T) / pi); HR(nii,:) = hr; % debilni definice sinc ! 
end
% priklad jednoho
subplot(413);
plot (zt,HR(14,zind),'b'); axis tight; hold on; stem(nT(14),xn(14)); hold off
% a ted vse, plus suma: 
xr= sum(HR); 
subplot(414);
plot (zt,HR(:,zind),'b',zt,xr(zind),'r'); axis tight; hold on; 
stem(nT(zindn),xn(zindn)); hold off
print -depsc FIG/rekoobd.eps

%%%%%%%%%%%%%%%%%%%%%%%
n = 0:100; xn = 5 * cos(0.025 * pi * n); stem(n,xn) 
print -depsc FIG/cos025pin.eps

n = 0:100; xn = 5 * cos(2.025 * pi * n); stem(n,xn) 
