clear all;
close all;

N = 5;
dt = 1e-1;
t = -N:dt:N-dt;
w = 2*pi*(0:1/length(t):1-1/length(t));

rect = @(x) double(abs(x)<1/2);

f = rect(t/3);
figure(1);plot(t,f);xlabel('temps');

F = fft(f);
figure(2);stem(abs(F));title('DFT coeffs');ylabel('module of the coefficients');

pause;

% approximation of the square window by fourier coeff (Gibbs phenomenon)
order = 15;

G = zeros(size(f));
G(1:ceil(order/2)) = F(1:ceil(order/2));
G(end-floor(order/2)+1:end) = F(end-floor(order/2)+1:end);

g = ifft(G);
figure(1);hold on;plot(t,real(g),'r-');hold off;
figure(2);hold on;stem(abs(G),'r-');hold off;
pause;

% phase lin?aire
module_H = ones(1,length(g));
angle_H = -2*w/dt;

figure(3);
subplot(2,1,1);plot(w,module_H,'k-');ylabel('module');xlabel('frequence');
subplot(2,1,2);plot(w,angle_H,'k-');ylabel('phase');xlabel('frequence');

H = module_H.*exp(1j*angle_H);

I = H.*G;

i = ifft(I);
figure(1);hold on;plot(t,real(i),'k-');xlabel('temps');ylabel('intensite');hold off;

pause;

% phase non-lin?aire
for distortion = 0.05:0.05:0.2
    module_H = ones(1,length(g));
    angle_H = -2*(w).^(1+distortion)/dt;
    
    figure(3);
    subplot(2,1,1);hold on;plot(w,module_H,'k--','LineWidth',2);hold off;
    subplot(2,1,2);hold on;plot(w,angle_H,'k--','LineWidth',2);hold off;
    
    H = module_H.*exp(1j*angle_H);
    
    I = H.*G;
    
    i = ifft(I);
    figure(1);hold on;plot(t,real(i),'k--','LineWidth',2);drawnow;hold off;
    
    pause;
end