%% metroda IIR
%Petr Procházka
close all

%% výběr signálu, svodu                
svod=1;                          %výběr svodu

y=xT(svod,:); 
fvz2 = 500;                     % vzorkovací frekvence
fz2 = 0;                        % stredni frekvence nepropustneho pasma
fi2=2*pi*fz2/fvz2;              % uhel
r2 = 0.985;                     % vzdalenost polu od pocatku
%%
% zadani koeficientu filtru
b2 = [1 -1];             % koeficienty citatele prenosove funkce (vypoctene z odvozeni z kruznice)
a2 = [1 -r2];            % koeficienty jmenovatele prenosove funkce
% vypocet frekvencni charakteristiky
N2 = 1000;
H2 = freqz(b2,a2,N2,'whole');           % vypocet 1000 bodu cele periody frekv.char.
b2=b2/max(H2);                          % timhle udelam jednotkovy prenos filtru    
H2 = freqz(b2,a2,N2,'whole');           % ale znovu musim nacist prenos s vynormovanymy koeficienty
% vypocet frekvencni osy
krok2 = fvz2/N2;
fosa2 = linspace(0,fvz2-krok2,N2);

% Pruchod signalu systemem
y2 = filtfilt( b2, a2, y );

%rozdílový signál
z=(xZ(svod,:)-y2);

cas= 0.002:0.002:10;            % Prepocet casove osy

%% vypočet SNR
orig=xZ(svod,:);
VykonSig=sum(orig.^2);
VykonSum=sum((y2-orig).^2);
SNR = 10*log( VykonSig/VykonSum )

%% vykreslime

figure
zplane(b2,a2); % rozlozeni nulovych bodu a polu
title('Rozložení nulových bodů a polů');
ylabel('Imaginární část')
xlabel('Reálná část')


figure
subplot(311)
plot(cas, xZ(svod,:))
title('originální signál')
xlabel('Čas [s]')
ylabel('napětí [uV]')
subplot(312);
plot( cas, y );
title('původní signál s navázaným trendem')
xlabel('Čas [s]')
ylabel('napětí [uV]')
subplot(313);
plot(cas, y2,'r');
title('vyfiltrovaný signál')
xlabel('Čas [s]')
ylabel('napětí [uV]')

figure
subplot(211)
plot(cas, xZ(svod,:))
hold on
plot(cas, y2,'r');
title('originální signál + vyfiltrovaný signál')
xlabel('Čas [s]')
ylabel('napětí [uV]')
subplot(212)
plot(cas, z)
title('rozdílový signál')
xlabel('Čas [s]')
ylabel('napětí [uV]')

spekt_puv=abs(fft(y));
spekt_vyf=abs(fft(y2));
spekt_orig=abs(fft(xZ(svod,:)));

delka_vyrezu_spektra=50;

figure
subplot(131);
stem((1:delka_vyrezu_spektra)/10,spekt_puv(1:delka_vyrezu_spektra));
title('              výřez spektra původního signálu')
axis([0 5 0 max(spekt_puv(1:delka_vyrezu_spektra))+10000])
xlabel('frekvence[Hz]')
ylabel('amplituda[uV/Hz]')
subplot(132);
stem((1:delka_vyrezu_spektra)/10,spekt_puv(1:delka_vyrezu_spektra));
title('              výřez spektra původního signálu')
axis([0 5 0 max(spekt_vyf(1:delka_vyrezu_spektra))+10000])
xlabel('frekvence[Hz]')
ylabel('amplituda[uV/Hz]')
subplot(133);
stem((1:delka_vyrezu_spektra)/10,spekt_vyf(1:delka_vyrezu_spektra),'r');
title('              výřez spektra vyfiltrovaného signálu')
axis([0 5 0 max(spekt_vyf(1:delka_vyrezu_spektra))+10000])
xlabel('frekvence[Hz]')
ylabel('amplituda[uV/Hz]')

% vykreslení spekter
delka_vyrezu_spektra1=100;

figure
subplot(211)
stem((1:delka_vyrezu_spektra1)/10,abs(spekt_puv(1:delka_vyrezu_spektra1)))
hold on
stem((1:delka_vyrezu_spektra1)/10,spekt_vyf(1:delka_vyrezu_spektra1),'r')
title('výřez spektra vyfiltrovaného a původního signálu s trendem')
xlabel('frekvence[Hz]')
ylabel('amplituda[uV/Hz]')
subplot(212)
stem((1:delka_vyrezu_spektra1)/10,spekt_orig(1:delka_vyrezu_spektra1))
hold on
stem((1:delka_vyrezu_spektra1)/10,spekt_vyf(1:delka_vyrezu_spektra1),'r')
title('výřez spektra originálního a vyfiltrovaného signálu')
xlabel('frekvence[Hz]')
ylabel('amplituda[uV/Hz]')

