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('----------------------------------------------------------------');
    disp(sprintf(' Fourierova analyza pro funkci zadanou tabulkou'));
    disp(' ');
    disp(sprintf(' |    x(k)    |    f(k)    |'));
    disp(sprintf(' ==========================='));
    for i=1:n
      disp(sprintf(' | %9.4f  | %9.4f  |',x(i),f(i)));
    end;  
    disp(sprintf(' | %9.4f  | %9.4f  |',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
	  disp(sprintf(...
	  ' koeficient A(%d) = %f    u bazove funkce phi(%d) = cos(2*pi*%d*x/%f-%f)',...
	  pocet/2,A(pocet/2),pocet-1,pocet/2,perioda,x(1)));
        else
	  disp(sprintf(...
	  ' koeficient B(%d) = %f    u bazove funkce phi(%d) = sin(2*pi*%d*x/%f-%f)',...
	  (pocet-1)/2,B((pocet-1)/2),pocet-1,(pocet-1)/2,perioda,x(1)));
	end;
      else
        disp(' ');
        disp(sprintf(' koeficient A(0) = %f    u bazove funkce phi(0) = 1',A0));
      end;
    end;	
	disp(' ');
	disp(' Aproximace je dana predpisem :');
	disp(' ');
	disp('----------------------------------------------------------------------');
        disp(sprintf('                 L'));
        disp(sprintf(...
	'  phi = A(0) + suma [ A(i)*phi(2i-1) + B(i)*phi(2i) ]'));
        disp(sprintf('                i=1'));
        disp(sprintf([...
	    '                                  pro pocet bazovych funkci N=2L+1']));
	disp('----------------------------------------------------------------------');
        disp(' ');
	disp('----------------------------------------------------------------------');
        disp(sprintf('                L-1'));
        disp(sprintf(...
	'  phi = A(0) + suma [ A(i)*phi(2i-1) + B(i)*phi(2i) ] + A(L)*phi(2L-1)'));
        disp(sprintf('                i=1'));
        disp(sprintf([...
	    '                                  pro pocet bazovych funkci N=2L']));
	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');
      popis1=sprintf(...
      'Aproximace trigonometrickym polynomem s %d bazovymi funkcemi  ',pocet);
      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
	  disp(sprintf(...
	  ' koeficient A(%d) = %f    u bazove funkce phi(%d) = cos(2*pi*%d*x/%f-%f)',...
	  pocet/2,A(pocet/2),pocet-1,pocet/2,perioda,x(1)));
        else
	  disp(sprintf(...
	  ' koeficient B(%d) = %f    u bazove funkce phi(%d) = sin(2*pi*%d*x/%f-%f)',...
	  (pocet-1)/2,B((pocet-1)/2),pocet-1,(pocet-1)/2,perioda,x(1)));
	end;
      else
        disp(' ');
        disp(sprintf(' koeficient A(0) = %f    u bazove funkce phi(0) = 1',A0));
      end;

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








