%PROPTF illustre les proprietes de la TF discrete
% voir : http://jmmeost.free.fr/Sismique/isis.html
% fig1 = http://jmmeost.free.fr/Sismique/proptf1.gif
% fig2 = http://jmmeost.free.fr/Sismique/proptf2.gif
% fig3 = http://jmmeost.free.fr/Sismique/proptf3.gif

clear all

nx = 4; dx = .25; LX = nx*dx; xx = 0:dx:(nx-1)*dx; %echantillonnage en x
dkx = 2*pi/LX; kxn = pi/dx; kx = 0:dkx:kxn; fx = kx/2/pi; fxn = kxn/2/pi; dfx = dkx/2/pi; %echantillonnage en kx et fx 
dxx = dx/10; xxx = 0:dxx:nx*dx; dkxx = dkx/10; kxx = -kxn:dkxx:kxn; fxx = kxx/2/pi; % representation des fonctions continues

for ifig = 1:2, figure(ifig), clf

for ik = 0:nx % nx+1 harmoniques, la derniere aliasee
    if ifig == 1, 
        P = exp(i*(ik-1)*dkx*xx); PX = exp(i*(ik-1)*dkx*xxx); end % harmoniques complexes
    if ifig ==2
        if ik<nx, P = cos(ik*dkx*xx); PX = cos(ik*dkx*xxx); else P = sin((ik-nx+1)*dkx*xx); PX = sin((ik-nx+1)*dkx*xxx); end  % harmoniques reelles
end

subplot(nx+1,4,1+4*ik), hold on
stem(xx,real(P),'r'); plot(xxx, real(PX),'r'); axis([0 LX -1 1]);
if (ik == 0) title('REEL( P )'); end, if (ik == nx) xlabel('X (M)'); end

subplot(nx+1,4,2+4*ik), hold on
stem(xx,imag(P),'r');  plot(xxx, imag(PX),'r'); axis([0 LX -1 1]);
if (ik == 0) title('IMAG( P )'); end, if (ik == 4) xlabel('X (M)'); end

P = fft(P); P = P/max(abs(P)); %TF et normalisation

subplot(nx+1,4,3+4*ik), hold on
stem(-fliplr(fx(2:nx/2)),real(P(nx/2+2:nx))); % -kxn+dkx <= kx <= -dkx (nx/2-1 valeurs)
stem(fx,real(P(1:nx/2+1))); % 0 <= kx <= kxn (nx/2+1 valeurs)
axis([-fxn fxn -1 1]);
if (ik == 0) title('REEL( TF(P) )'); end, if (ik == nx) xlabel('FX (CYCLE/M)'); legend('ALIASING'); end

subplot(nx+1,4,4+4*ik), hold on
stem(-fliplr(fx(2:nx/2)),imag(P(nx/2+2:nx))); stem(fx,imag(P(1:nx/2+1)));
axis([-fxn fxn -1 1]);
if (ik == 0) title('IMAG( TF(P) )'); end, if (ik == nx) xlabel('FX (CYCLE/M)'); end

end
end

figure(3), clf

for ip = 0:nx-1
    
P = zeros(1,nx); P(ip+1) =1; % diracs reels

subplot(nx,4,1+4*ip), hold on
stem(xx,real(P),'r'); axis([0 LX -1 1]);
if (ip == 0) title('REEL(P)'); end, if (ip == nx-1) xlabel('X (M)'); end

subplot(nx,4,2+4*ip), hold on
stem(xx,imag(P),'r');  axis([0 LX -1 1]);
if (ip == 0) title('IMAG( P )'); end, if (ip == nx-1) xlabel('X (M)'); end

P = fft(P); P = P/max(abs(P)); %TF et normalisation
PX = exp(-i*ip*dx*kxx); % TF continue

subplot(nx,4,3+4*ip), hold on
stem(-fliplr(fx(2:nx/2)),real(P(nx/2+2:nx))); stem(fx,real(P(1:nx/2+1))); 
plot(fxx,real(PX),'b'); 
plot(-fx(2),real(P(2)),'r*',fx(2),real(P(2)),'r*'); % Reel (TF(P reel)) symetrique
axis([-fxn fxn -1 1]);
if (ip == 0) title('REEL( TF(P) )'); end, if (ip == nx-1) xlabel('FX (CYCLE/M)'); end

subplot(nx,4,4+4*ip), hold on
stem(-fliplr(fx(2:nx/2)),imag(P(nx/2+2:nx))); stem(fx,imag(P(1:nx/2+1))); 
plot(fxx,imag(PX),'b'); axis([-fxn fxn -1 1]);
plot(-fx(2),-imag(P(2)),'r*',fx(2),imag(P(2)),'r*'); % Imag (TF(P reel)) antisymetrique
if (ip == 0) title('IMAG( TF(P) )'); end, if (ip == nx-1) xlabel('FX (CYCLE/M)'); end

end