function fourierova_analyza(zobraz_vse,x,f,stupen)
%
% provede diskretni Fourierovu analyzu 
% pro periodickou funkci zadanou tabulkou 
%
% stupen ... pocet bazovych funkci
%            ( pokud neni zadan stupen vezme 
%              maximalni mozny pocet bazovych funkci
%              odpovidajici poctu zadanych bodu )     
% uvazujeme ekvidistantni uzly x, 
% pricemz posledni bod nezadavame 
% (funkcni hodnota je stejna jako u prvniho bodu
%  a x-ovou souradnici lehko urcime, nebot 
%  uvazujeme ekvidistantni rozlozeni uzlu)

clc;
epsilon=1e-8;
pocet_deleni=500;
pom=diff(x);
pom=pom/pom(1);
test=0;
for i=1:max(size(pom))
  if abs(pom(i)-1)>epsilon
    test=1;
  end
end
n=max(size(x));
deleni=(x(n)-x(1))/pocet_deleni;


if exist('stupen')
  if stupen>n
    disp(' ');
    disp(' Pozor! Pocet bazovych funkci nesmi byt vetsi nez pocet zadanych bodu');
    disp(' ');
    stupen=n;
  end
else
  stupen=n;
end  


if test==1
  disp(' ');
  disp(' Pozor! Zadany vektor x nema ekvidistantni uzly. ');
  disp(' ');
else  
    perioda=n*(x(n)-x(1))/(n-1);
    disp('----------------------------------------------------------------');
    fprintf(' Fourierova analyza pro funkci zadanou tabulkou\n');
    disp(' ');
    fprintf(' |    x(k)    |    f(k)    |\n');
    fprintf(' ===========================\n');
    for i=1:n
      fprintf(' | %9.4f  | %9.4f  |\n',x(i),f(i));
    end 
    fprintf(' | %9.4f  | %9.4f  |\n',x(1)+perioda,f(1));
    disp(' ');
    C=fft(f);    
    max_stupen=round((length(C)-1)/2);
    A0=C(1)/n;
    if abs(A0)<epsilon
      A0=0;
    end  
    N=min([round((stupen-1)/2),max_stupen]);
    for k=1:N
      A(k) = 2*real(C(k+1))/n;
      if abs(A(k))<epsilon
        A(k)=0;
      end  
      B(k) = -2*imag(C(k+1))/n;
      if abs(B(k))<epsilon
        B(k)=0;
      end  
    end
    if C(N+1)~=0 && ~mod(length(C),2)
      A(N)=A(N)/2;
    end    
    disp(' Stiskni klavesu'); 
    pause;
    figure(1);
    set(figure(1),'Units','centimeters');
    set(figure(1),'Position',[22 6 18 15]);
    clf;
    hold on;
    zoom on;
 if zobraz_vse~=1
    xx=x(1):deleni:x(1)+perioda;
    l=0;
    fmin=min(f);
    fmax=max(f);
    for s=0:deleni:perioda
      l=l+1;
      ff(l)=A0;
      for m=1:N
	    ff(l)=ff(l)+A(m)*cos(2*pi/perioda*m*s);
	    if 2*m<=stupen-1 
	      ff(l)=ff(l)+B(m)*sin(2*pi/perioda*m*s);
        end
      end	
    end
    fmin=min([fmin, min(ff)]);
    fmax=max([fmax, max(ff)]);
    axis([x(1)-perioda/10 x(1)+perioda*11/10 fmin-(fmax-fmin)/10 fmax+(fmax-fmin)/10]);       
    plot(xx,ff,'b');
    plot([x,x(1)+perioda],[f,f(1)],'ro');
    for pocet=1:min([stupen, max_stupen*2+1])
      z=mod(pocet,2);
      if pocet>1
	    if ~z
	      fprintf(...
	  ' koeficient A(%d) = %f    u bazove funkce phi(%d) = cos(2*pi*%d*x/%f-%f)\n',...
	  pocet/2,A(pocet/2),pocet-1,pocet/2,perioda,x(1));
        else
	  fprintf(...
	  ' koeficient B(%d) = %f    u bazove funkce phi(%d) = sin(2*pi*%d*x/%f-%f)\n',...
	  (pocet-1)/2,B((pocet-1)/2),pocet-1,(pocet-1)/2,perioda,x(1));
        end
      else
        disp(' ');
        fprintf(' koeficient A(0) = %f    u bazove funkce phi(0) = 1\n',A0);
      end
    end	
	disp(' ');
	disp(' Aproximace je dana predpisem :');
	disp(' ');
	disp('----------------------------------------------------------------------');
        fprintf('                 L\n');
        fprintf(...
	'  phi = A(0) + suma [ A(i)*phi(2i-1) + B(i)*phi(2i) ]\n');
        fprintf('                i=1\n');
        fprintf(...
	    '                                  pro pocet bazovych funkci N=2L+1\n');
	disp('----------------------------------------------------------------------');
        disp(' ');
	disp('----------------------------------------------------------------------');
        fprintf('                L-1\n');
        fprintf(...
	'  phi = A(0) + suma [ A(i)*phi(2i-1) + B(i)*phi(2i) ] + A(L)*phi(2L-1)\n');
        fprintf('                i=1\n');
        fprintf(...
	    '                                  pro pocet bazovych funkci N=2L\n');
	disp('----------------------------------------------------------------------');
  else  
    stupen=min([2*max_stupen+1,stupen]);
    xx=x(1):deleni:x(1)+perioda;
    fmin=min(f);
    fmax=max(f);
    for pocet=1:stupen 
      l=0;
      for s=x(1):deleni:x(1)+perioda
        l=l+1;
        ff(l)=A0;
        for m=1:floor(pocet/2) 
	  ff(l)=ff(l)+A(m)*cos(2*pi/perioda*m*s);
          if 2*m<=stupen-1  && 2*m<pocet
	    ff(l)=ff(l)+B(m)*sin(2*pi/perioda*m*s);
    	  end
        end
      end
      fmin=min([fmin, min(ff)]);
      fmax=max([fmax, max(ff)]);
    end
    for pocet=1:stupen 
      l=0;
      for s=0:deleni:perioda
        l=l+1;
        ff(l)=A0;
        for m=1:floor(pocet/2) 
	  ff(l)=ff(l)+A(m)*cos(2*pi/perioda*m*s);
          if 2*m<=stupen-1  && 2*m<pocet
	    ff(l)=ff(l)+B(m)*sin(2*pi/perioda*m*s);
    	  end
        end	
      end
      axis([x(1)-perioda/10 x(1)+perioda*11/10 fmin-(fmax-fmin)/10 fmax+(fmax-fmin)/4]);       
      plot(xx,ff,'b');
      plot([x,x(1)+perioda],[f,f(1)],'ro');
      if pocet==1
        popis1=sprintf(...
       'Aproximace trigonometrickym polynomem s 1 bazovou funkci      ',pocet);
      else
       popis1=sprintf(...
       'Aproximace trigonometrickym polynomem s %d bazovymi funkcemi  ',pocet);
      end    
      popis2=sprintf(...
      'aproximovane body                                              ');
      texty(1,:)=popis1(1:60); 
      texty(2,:)=popis2(1:60); 
      legend(texty);

      z=mod(pocet,2);
      if pocet>1
	if ~z
	  fprintf(...
	  ' koeficient A(%d) = %f    u bazove funkce phi(%d) = cos(2*pi*%d*x/%f-%f)\n',...
	  pocet/2,A(pocet/2),pocet-1,pocet/2,perioda,x(1));
        else
	  fprintf(...
	  ' koeficient B(%d) = %f    u bazove funkce phi(%d) = sin(2*pi*%d*x/%f-%f)\n',...
	  (pocet-1)/2,B((pocet-1)/2),pocet-1,(pocet-1)/2,perioda,x(1));
    end
      else
        disp(' ');
        fprintf(' koeficient A(0) = %f    u bazove funkce phi(0) = 1\n',A0);
      end

      if pocet<stupen
        pause;  
        clf;
        hold on;
      else
	disp(' ');
	disp(' Aproximace je dana predpisem :');
	disp(' ');
	disp('----------------------------------------------------------------------');
        fprintf('                 L\n');
        fprintf(...
	'  phi = A(0) + suma [ A(i)*phi(2i-1) + B(i)*phi(2i) ]\n');
        fprintf('                i=1\n');
        fprintf(...
	    '                                  pro pocet bazovych funkci N=2L+1\n');
	disp('----------------------------------------------------------------------');
        disp(' ');
	disp('----------------------------------------------------------------------');
        fprintf('                L-1\n');
        fprintf(...
	'  phi = A(0) + suma [ A(i)*phi(2i-1) + B(i)*phi(2i) ] + A(L)*phi(2L-1)\n');
        fprintf('                i=1\n');
        fprintf(...
	    '                                  pro pocet bazovych funkci N=2L\n');
	disp('----------------------------------------------------------------------');
      end	
    end
 end  
end  








