% Ilustrace algoritmu myfft a myifft

% Ukazkove signaly
n=8; N=2^n;
%t=(fix(-N/2):fix(N/2-1));
t=(0:N-1)';
w1=2*pi/N; w2=2*pi/N*16; w3=2*pi/N*4;
f1=sin(w1*t)+sin(w2*t)+sin(w3*t);
f2=sin(w1*t);
f2(fix(N/4):fix(3*N/4))=f2(fix(N/4):fix(3*N/4))+0.3*sin(w2*t(fix(N/4):fix(3*N/4)));
f2(fix(N/2):N-1)=f2(fix(N/2):N-1)+0.3*sin(w3*t(fix(N/2):N-1));

% Prima transformace
F1=fft(f1);
myF1=myfft(f1);
F2=fft(f2);
myF2=myfft(f2);


% Vynulovani "nulovych" koeficientu (nenuly v dusledku zaokrouhlovani)
Eps=1e-10; % Cisla mensi nez Eps povazujeme za nulove

Re=real(F1); iRe=abs(Re)>=Eps;
Im=imag(F1); iIm=abs(Im)>=Eps;
F1=Re.*iRe+j*Im.*iIm;

Re=real(myF1); iRe=abs(Re)>=Eps;
Im=imag(myF1); iIm=abs(Im)>=Eps;
myF1=Re.*iRe+j*Im.*iIm;

Re=real(F2); iRe=abs(Re)>=Eps;
Im=imag(F2); iIm=abs(Im)>=Eps;
F2=Re.*iRe+j*Im.*iIm;

Re=real(myF2); iRe=abs(Re)>=Eps;
Im=imag(myF2); iIm=abs(Im)>=Eps;
myF2=Re.*iRe+j*Im.*iIm;


%vypocet ampl. a faz. spektra
F1_am=abs(F1)/(N/2);
myF1_am=abs(myF1)/(N/2);
F2_am=abs(F2)/(N/2);
myF2_am=abs(myF2)/(N/2);

F1_ph=angle(F1);
myF1_ph=angle(myF1);
F2_ph=angle(F2);
myF2_ph=angle(myF2);

% Zpetna transformace pomoci myifft
myf1_=real(myifft(myF1)); % Eliminace komplexni casti pomoci funkce real
myf2_=real(myifft(myF2));

% Graficke zobrazeni spekter
figure(1);
  subplot(2,2,1);
    hold on
    axis([0 N/2-1 min(F1_am(1:N/2)) max(F1_am(1:N/2))]);
    title('Amplitudové spektrum f_1 pomocí fft','FontName','@MS Mincho');
    bar(0:N/2-1,F1_am(1:N/2));
  subplot(2,2,2);
    hold on
    axis([0 N/2-1 min(myF1_am(1:N/2)) max(myF1_am(1:N/2))]);
    title('Amplitudové spektrum f_1 pomocí myfft','FontName','@MS Mincho');
    bar(0:N/2-1,myF1_am(1:N/2));
  subplot(2,2,3);
    hold on
    axis([0 N/2-1 min(F1_ph(1:N/2)) max(F1_ph(1:N/2))]);
    title('Fázové spektrum f_1 pomocí fft','FontName','@MS Mincho');
    bar(0:N/2-1,F1_ph(1:N/2));
  subplot(2,2,4);
    hold on
    axis([0 N/2-1 min(myF1_ph(1:N/2)) max(myF1_ph(1:N/2))]);
    title('Fázové spektrum f_1 pomocí myfft','FontName','@MS Mincho');
    bar(0:N/2-1,myF1_ph(1:N/2));
  pause;
  clf;
  subplot(2,2,1);
    hold on
    axis([0 N/2-1 min(F2_am(1:N/2)) max(F2_am(1:N/2))]);
    title('Amplitudové spektrum f_2 pomocí fft','FontName','@MS Mincho');
    bar(0:N/2-1,F2_am(1:N/2));
  subplot(2,2,2);
    hold on
    axis([0 N/2-1 min(myF2_am(1:N/2)) max(myF2_am(1:N/2))]);
    title('Amplitudové spektrum f_2 pomocí myfft','FontName','@MS Mincho');
    bar(0:N/2-1,myF2_am(1:N/2));
  subplot(2,2,3);
    hold on
    axis([0 N/2-1 min(F2_ph(1:N/2)) max(F2_ph(1:N/2))]);
    title('Fázové spektrum f_2 pomocí fft','FontName','@MS Mincho');
    bar(0:N/2-1,F2_ph(1:N/2));
  subplot(2,2,4);
    hold on
    axis([0 N/2-1 min(myF2_ph(1:N/2)) max(myF2_ph(1:N/2))]);
    title('Fázové spektrum f_2 pomocí myfft','FontName','@MS Mincho');
    bar(0:N/2-1,myF2_ph(1:N/2));
  pause;
  clf;
  subplot(2,2,1);
    hold on
    axis([min(t) max(t) min(f1) max(f1)]);
    title('Pùvodní signál f_1','FontName','@MS Mincho');
    plot(t,f1,'b-');
  subplot(2,2,2);
    hold on
    axis([min(t) max(t) min(f2) max(f2)]);
    title('Pùvodní signál f_2','FontName','@MS Mincho');
    plot(t,f2,'b-');
  subplot(2,2,3);
    hold on
    axis([min(t) max(t) min(myf1_) max(myf1_)]);
    title('Zrekonstruovaný signál f_1^* pomoci myifft','FontName','@MS Mincho');
    plot(t,myf1_,'g-');
  subplot(2,2,4);
    hold on
    axis([min(t) max(t) min(myf2_) max(myf2_)]);
    title('Zrekonstruovaný signál f_2^* pomoci myifft','FontName','@MS Mincho');
    plot(t,myf2_,'g-');
  pause;
close(1);
