function DiscreteFourierTransform

i=20;

N=i*100;
Resolution = 1e-4;
TimeDomain = i*pi;
SaveName = strcat('DFT_nonlin',num2str(i),'pi');


GlobalError=zeros(N,1);
for k=N:N
    
    ti = 0:TimeDomain/(2*k+1):2*k*TimeDomain/(2*k+1);
    fi = zeros(size(ti));
    for i=1:length(ti)
        fi(i) = f(ti(i));
    end
    
    [freqs, amps, phases, coefs] = real_fft( ti, fi );
    
    MeshX = 0:Resolution:TimeDomain;
    FV = zeros(size(MeshX));
    for i=1:length(MeshX)
        FV(i) = f(MeshX(i));
    end
    
    FT = coefs(1) * exp(1i * 0 * 2*pi / TimeDomain * MeshX) / (2*k+1);
    for i=1:k
        FT = FT + coefs(i+1)     * exp(1i *  i * 2*pi / TimeDomain * MeshX) / (2*k+1);
        FT = FT + coefs(end-i+1) * exp(1i * -i * 2*pi / TimeDomain * MeshX) / (2*k+1);
    end
    
    fig = figure(1);
        set(fig,'Position',[680 560 1120 420]);
    
    subplot(1,2,1); hold on;
        set(gca,'NextPlot','replace');
            pFV = plot(MeshX,FV);
        set(gca,'NextPlot','add');
            pFT = plot(MeshX,FT);
            pTI = plot(ti,fi);
            set(gca,'Box','on','XGrid','on','YGrid','on','XLim',[0 TimeDomain]);
            set(gca,'FontSize',12,'FontName','Times','FontWeight','bold');
        xlabel('t(s)','FontSize',16,'FontName','Times');
        ylabel('f(t)','FontSize',16,'FontName','Times');
        title(strcat('N=',num2str(k)),'FontSize',20,'FontName','Times');
            set(pFV,'LineWidth',2,'Color',[0 0 0]);
            set(pFT,'LineWidth',2,'Color',[1 0 0]);
            set(pTI,'LineStyle','none','Color',[0 0 0],'Marker','.','MarkerSize',10); hold off;
    subplot(1,2,2);
        spectra = bar(2*pi*freqs,amps);
            set(gca,'Box','on','XGrid','on','YGrid','on','XLim',[0 12],'YLim',[0 5]);
            set(gca,'FontSize',12,'FontName','Times','FontWeight','bold');
        xlabel('\omega_i(rad/s)','FontSize',16,'FontName','Times');
        ylabel('A(\omega)','FontSize',16,'FontName','Times');
            set(spectra,'BarWidth',0.4,'FaceColor',[0 0 0]);
    pause(0.1)
    saveas(fig, strcat(SaveName,num2str(k)), 'png');
    
    ErrorFunction    = abs( FV - FT );
    GlobalError(k+1) = max(ErrorFunction);
    
%     pause
end
        
% ListOfN=0:1:N;
% 
% figure(3);
% err=plot(ListOfN,GlobalError);
%     set(gca,'Box','on','XGrid','on','YGrid','on','XLim',[0 N],'YScale','log');
%     set(gca,'FontSize',12,'FontName','Times','FontWeight','bold');
%         xlabel('N','FontSize',16,'FontName','Times');
%         ylabel('E','FontSize',16,'FontName','Times');
%     set(err,'LineWidth',2,'Color',[0 0 0]);


% saveas(fig, strcat(SaveName,num2str(N)), 'png');






function y=f(x)

% y = 3./(5-4*cos(x));
% y = abs(-x + pi);
% y = sign(x - pi);
% y = 3*cos(x) + 5*sin(5*x);
% y = sin(9*x);

% if x<2*pi
%     y = sin(x);
% else
%     y = sin(5*x);
% end

% y = (1+0.2*sin(0.2*x+2)) * sin(3*x);
% y = sin( (3+0.2*sin(0.2*x)) * x);
% y = sin( (3+0.01*sin(0.2*x)) * x);

% y = sin(0.3*x) + 1.3*sin(3*x) + 0.36*sin(2*x) + 2.1*sin(7*x) + 3.7*sin(10*x) + 1*sin(0.6*x) + 0.3*sin(4*x);

y = sin(0.3*x) + exp(1.3*sin(3*x)) + 0.36*sin(2*x) + 2.1*sin(7*x) + 3.7*sin(10*x) + 1*sin(0.6*x) + 0.3*sin(4*x);

% y = exp( sin(3*x) );

function y=Cos0(x)
y=1 + 0*x;

function y=CosN(x,N)
y=cos(N*x);

function y=SinN(x,N)
y=sin(N*x);


function [freqs, amps, phases, coefs] = real_fft( t, signal )

timestep = t(2)-t(1);
N = length(signal);
N_half = ceil((N+1)/2);
freqs(1:N_half,1) = 1/timestep * (0:N_half-1)/N; % 1/timestep = mintavételezési frekvencia

coefs = fft(signal);

amps(1:N_half,1) = 2/N*abs( coefs(1:N_half) );
amps(1) = amps(1)/2;
if (mod(N,2)==0) amps(N_half)=amps(N_half)/2; end;
phases(1:N_half,1) = angle(coefs(1:N_half));