%=================================================================
%         FILTRACE METODOU POSTUPNE POLYNOMIALNI REGRESE
% - vypocet vahovych koeficentu (impulzni odezvy linearniho filtru)
% - vlastni vyhlazeni provedenim konvoluce (klouzavy prumer)
%=================================================================

% (C) 1996-98, V. Vesely, Masaryk University of Brno (Czech Republic)

format compact;
format long
echo on
% Vypocet vahovych koeficientu filtru delky m=2*q+1 metodou 
% lokalni regrese polynomem stupne d:
format rational
  d=2; p=d+1; m=5; q=fix(m/2); D=vander(-q:q); D=D(:,m-p+1:m); H=inv(D'*D)*D';
  D              % Matice planu pro lokalni regresi (t=0 = stred vahoveho okna) 
  for j=1:q
    w=D(j,:)*H % vahove okno pro levy usek dat k=j=1,...,q   
  end
pause;
%---------------------------------------------------------------
w=H(p,:)     % vahove okno pro stredni usek dat k=q+1,...,n-q (n=delka dat)
w=D(q+1,:)*H % je totez, nebot prostredni radek D je nulovy az na 1 v posl. sl.
pause;
%---------------------------------------------------------------
  for j=q+2:m
    w=D(j,:)*H % vahove okno pro pravy usek dat k=n-m+j, tj. k=n-q+1,...,n   
  end
pause;
%---------------------------------------------------------------
% ZOPAKUJME TOTEZ PRO o 1 vyssi stupen polynomu d=3;
  d=3; p=d+1; m=5; q=fix(m/2); D=vander(-q:q); D=D(:,m-p+1:m); H=inv(D'*D)*D';
  for j=1:q
    w=D(j,:)*H % vahove okno pro levy usek dat k=j=1,...,q   
  end
pause;
%---------------------------------------------------------------
w=H(p,:)     % vahove okno pro stredni usek dat k=q+1,...,n-q (n=delka dat)
% VSIMNETE SI, ZE VAHA PRO STREDNI USEK JE STEJNA JAKO pro d=2;
% Da se ukazat, ze toto plati obecne pro kazde stupne d a d+1, kde d je sude.
% Pro okrajove vahy to vsak jiz neplati. 
% Vsimnete si take techto symetrii:
% - Stredove vahy jsou symetricke vzhledem ke stredu
% - Okrajove vahy pro j-ty bod a symetricky (n+1-j)-ty bod jsou stredove otocene.  
pause;
%---------------------------------------------------------------
  for j=q+2:m
    w=D(j,:)*H % vahove okno pro pravy usek dat k=n-m+j, tj. k=n-q+1,...,n   
  end
pause;
%---------------------------------------------------------------
% vaha pro predikci o r bodu do minulosti, napr. r=1:
  r=1;
  wl1=((-q-r).^(d:-1:0))*H
% vaha pro predikci o r bodu do budoucnosti:
  wp1=((q+r).^(d:-1:0))*H
% Opet lze pozorovat symetrii v otoceni
pause;
%===========================================================================
% PROVEDEME NYNI VLASTNI VYHLAZENI POMOCI KONVOLUCE
% GENEROVANI SIMULACNICH DAT
  n = 100;
  t = (0:n-1).';
  y = 2 + sin(2*pi*t./n);  % Sinusovy trend (1 perioda)
  sigma = 0.2;
  x = y + sigma*randn(n,1);
  figure; plot(t,x,'w'); hold on; plot(t,y,'r','LineWidth',2);
  title('SIMULACE SINUSOVKY KONTAMINOVANE BILYM SUMEM');
pause;
%---------------------------------------------------------------
% Z predchoziho postupu je jasne, ze vsechny vahy (bez predikce)
% ziskame jako radky teto matice:
  d=3; p=d+1; m=5; q=fix(m/2); D=vander(-q:q); D=D(:,m-p+1:m); H=inv(D'*D)*D';
  FIR=D*H
  yl = FIR(1:q,:)*x(1:m);        % Vyhlazeni leveho okraje
  ys = conv(FIR(q+1,:),x);
  ys = ys(m:n);                  % Vyhlazeni stredni casti
  yp= FIR(q+2:m,:)*x(n-m+1:n);   % Vyhlazeni praveho okraje
  yy = [yl;ys;yp];               % Celkove vyhlazeni
  plot(t,yy,'g','LineWidth',2);  % Prikreslime do grafu
  title(['SINUSOVKA FILTROVANA S OKNEM DELKY ',num2str(m)]);
  hold off;
pause;
%---------------------------------------------------------------
% Totez provedeme pohodlneji pomoci procedury convfilt
  yyy = convfilt(FIR,x);
  max(abs(yy-yyy))     % Vysledek je tyz
pause;
%---------------------------------------------------------------
% Vyhlazeni je zjevne nedostatecne.
% Zkusime zvetsit delku filtru na 21:
  d=3; p=d+1; m=21; q=fix(m/2); D=vander(-q:q); D=D(:,m-p+1:m); H=inv(D'*D)*D';
  FIR=D*H;
  yy = convfilt(FIR,x);
  figure; plot(t,x,'w'); hold on; plot(t,y,'r','LineWidth',2);
  plot(t,yy,'g','LineWidth',2);  % Prikreslime do grafu
  title(['SINUSOVKA FILTROVANA S OKNEM DELKY ',num2str(m)]);
  hold off;
pause;
%---------------------------------------------------------------
% Jeste zkusime zvetsit  delku filtru na 31:
  d=3; p=d+1; m=31; q=fix(m/2); D=vander(-q:q); D=D(:,m-p+1:m); H=inv(D'*D)*D';
  FIR=D*H;
  yy = convfilt(FIR,x);
  figure; plot(t,x,'w'); hold on; plot(t,y,'r','LineWidth',2);
  title('SIMULACE SINUSOVKY KONTAMINOVANE BILYM SUMEM');
  plot(t,yy,'g','LineWidth',2);  % Prikreslime do grafu
  title(['SINUSOVKA FILTROVANA S OKNEM DELKY ',num2str(m)]);
  hold off;
pause;
%---------------------------------------------------------------
% Nyni to jiz vypada dobre.
% ZKUSIME PREDIKCI
% a) DO MINULOSTI O 1 KROK:
  format long;
  wl1=((-q-1).^(d:-1:0))*H;
  yym1 = wl1*x(1:m)            % Predikovana hodnota v case t=-1
  ym1  = 2 + sin(2*pi*(-1)./n) % Presny sinusovy trend v case t=-1
% b) DO BUDOUCNOSTI O 1 KROK:
  wp1=((q+1).^(d:-1:0))*H;
  yyp1 = wp1*x(n-m+1:n)     % Predikovana hodnota v case t=n
  yp1  = 2 + sin(2*pi*n./n) % Presny sinusovy trend v case t=n
% Odhad je prijatelny. Dale do minulosti, resp. budoucnosti by se ovsem
% kvalita predikce zhorsovala.
%---------------------------------------------------------------------------
% PROSTUDUJTE ALGORITMUS PROCEDUR 'convfilt' pomoci prikazu 'type'.
% Zkuste vyhlazovani simulovanych dat s vyssim stupnem: d=4,6,8,...
%================================== KONEC ==================================
echo off
