% close all
% clear all
% clc

%load hdpcz;
function [yp2] = CHF(hdpcz)
%hdpcz=ea12;
outctv=hpf(hdpcz,1600);
rozctv=hdpcz-outctv;
furctv=abs(fft(rozctv));

osaf1=1:1:length(furctv)/2;       %vypocet frakci pi
osaf1=osaf1./length(osaf1);    %vypocet frakci pi
rpul = furctv(1:ceil(length(furctv)/2));

% figure(1);plot(rozctv);xlim([1, length(hdpcz)]);%grid on    %cyklicka slozka po filtraci HP
% figure(2);plot(osaf1,rpul,'black');%grid on   %spektrum pofiltraci HP
% hold on
%    % line([2/32 2/6],[max(rpul) max(rpul)])
%     line([2/32 2/32],[0 max(rpul)])
%     line([2/6 2/6],[0 max(rpul)])
%     %pause
% hold off
%---------------------------------------------------------
%filtr CH-F
x = hdpcz;  %vstupni posloupnost

T=length(x);

% vypocet koeficientu Bj, B0 az BT-1 ulozenych do vektoru B[B0,B1,B2,...]
a = 2*pi/32;
b = 2*pi/6;

B = zeros(1,T);
for j=1:T-1
    B(1,j+1)=(sin(j*b)-sin(j*a))/(j*pi);
end
B(1,1)=(b-a)/pi; 

%---------------------------------------------------------------
%vypocet prvni casti vyrazu, vystupem vektor yp1[0, 0, yp1(t), 0, 0]

yp1 = zeros(1,T);
BB=zeros(T,T);      %prenos

pom = 0;                            %pomocna promenna
for i=3:T-2                 % cyklus pro vsechny cleny yp1
    k = zeros(1,T);                 %pomocny vektor
    for t=3+pom:T-1             % vypocet jednoho clenu yp1 
        k(t) = x(t)*B(t-2-pom);
        BB(i,t)= B(t-2-pom);    %prenos
    end
    pom = pom+1;
    yp1(i) = sum(k);
end
yp1;
%---------------------------------------------------------------
%vypocet druhe a ctvrte casti 

%vypocet Bp,vystupem je vektor Bp[0, Bp(T-t), 0, 0, 0]

Bp = zeros(1,T);        
for t=3:T-2
    k = zeros(1,T);     %pomocny vektor
    Bs = 0;             %pomocna promenna
    for j=1:T-t-1
        k(j) = B(j+1);          %prislusne koeficienty "Bj"
    end
    Bs=sum(k);                  %suma prislusnych koeficientu "Bj"
    Bp(T-t) = -1/2*B(1)-Bs;     
end
%Bp

%vypocet druheho a ctvrteho clenu (jejich soucet)

yp24 = zeros(1,T);
for t=3:T-2
    yp24(t) = Bp(T-t)*x(T) + Bp(t-1)*x(1);
    BB(t,1)=Bp(t-1);     %prenos
    BB(t,T)=Bp(T-t);     %prenos
end
yp24;
%---------------------------------------------------------------
%vypocet treti casti

yp3 = zeros(1,T);

pom = 0;                            %pomocna promenna
for i=3:T-2                         % cyklus pro vsechny cleny yp1
    k = zeros(1,T);                 %pomocny vektor
    for t=3+pom:3+pom+pom                       % vypocet jednoho clenu yp3 
        k(t) = x(3+pom-1-(t-3-pom))*B(t-1-pom);
        BB(i,3+pom-1-(t-3-pom))=B(t-1-pom);      %prenos
    end
    pom = pom+1;
    yp3(i) = sum(k);
end

yp3;
%---------------------------------------------------------------
%celkovy vypocet

yp = zeros(1,T);

yp = yp1 + yp24 + yp3;

%---------------------------------------------------------------
%postprocessing
%casova oblast
% figure (33)
% plot(yp');
% xlim([1, length(yp)]);  %nastaveni rozsahu osy x[od,do]
%grid on

%kmitoctova oblast
r=abs(fft(yp));

osaf=1:1:ceil(length(r)/2);       %vypocet frakci pi
osaf=osaf./length(osaf);    %vypocet frakci pi
   
figure (44)
rpul = r(1:ceil(length(r)/2));
plot(osaf,rpul,'black');            %frakce pi na ose x
hold on
   % line([2/32 2/6],[max(rpul) max(rpul)])
    line([2/32 2/32],[0 max(rpul)])
    line([2/6 2/6],[0 max(rpul)])
    %pause
hold off
%plot(rpul);
%xlim([1, length(r)/2]);
%grid on;
%--------------------------------------------------------------
%vykresleni prenosu - pro jednotlive t odblokovat pause
% for t=3:T-2
%     t   
%     figure(55)
%     [HH,ww]=freqz(BB(t,:),1,0:pi/1000:pi);
%     plot(ww./pi,abs(HH),'black');
%     hold on
%     line([2/32 2/6],[1 1])
%     line([2/32 2/32],[0 1])
%     line([2/6 2/6],[0 1])
%     %pause
%     hold off
% end

yp2=yp(1,3:length(yp)-2);
plot(yp2)