%=================================================================
% FILTRACE METODOU POSTUPNE POLYNOMIALNI REGRESE POMOCI UNIVERZALNI
%                          PROCEDURY 'POLYFILT'
%=================================================================

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

format compact;
format long
echo on
% Funkce polyfilt umoznuje pohodlne reseni vseho, co souvisi s postupnou
% polynomialni regresi:
% 1) Provedeni operace klouzaveho prumeru na zadane siti vcetne vypoctu
%    konfidencnich intervalu pro vyhlazene i predikovane hodnoty
% 2) Vraceni operatorove matice FIR, pomoci niz muzeme realizovat
%    vlastni vyhlazeni pozdeji jako y=FIR*x, kde x je vektor puvodnich dat
%    a y vyhlazeny prubeh na zadane siti.
% 3) Je-li vystupni sit pro y totozna se siti pro x, pak radky FIR
%    reprezentuji vahy porade pro levy okraj, stredovy usek a pravy okraj tak
%    jak byly explicitne pocitany v crcv7.
%-----------------------------------------------------------------
% PROVEDTE NYNI help polyfilt A PROSTUDUJTE OVLADANI TETO PROCEDURY
%-----------------------------------------------------------------
keyboard;

%=================================================================
% VYGENERUJEME SIMULACNI DATA JAKO V crcv7:
  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;
format rational
%---------------------------------------------------------------
% Misto primeho vypoctu vah jako v crcv7, tj.:
  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';
  FIR1=D*H;
% provedeme jednoduse
  FIR2=polyfilt(d,m)
  max(abs(FIR1-FIR2)) % Vidime, ze matice vah jsou totozne
pause;
format long
%---------------------------------------------------------------
% Vlastni vyhlazeni pak provedeme bud pomoci 'convfilt' jako v crcv7, kde
% prvni parametr FIR vypocteme rovnou pomoci 'polyfilt':
  y1 = convfilt(polyfilt(d,m),x);
% Nebo ekvivalentne toto vsechno rovnou udela 'polyfilt' takto:
  y2 = polyfilt(d,m,x);
  max(abs(y1-y2))     %  Vidime opet shodne vysledky z obou postupu
pause;
%---------------------------------------------------------------
% Prime pouziti 'polyfilt' pro vyhlazeni vsak umoznuje ziskat navic
% i polosirky konfidencnich intervalu zadanim druheho vystupniho parametru E:
  [ys,E] = polyfilt(d,m,x);
% Chceme-li vsechny vystupni udaje a soucasne graficke znazorneni vysledku
% vyhlazeni, bude volani vypadat takto:
  [ys,E,r,s,bE] = polyfilt(d,m,x,t,[],[],figure);
pause;
%---------------------------------------------------------------
% V crcv7 jsme zjistili, ze vyhovuje delka m=31:
  m=31;
  figure; plot(t,y,'r','LineWidth',2); hold on
  [ys,E,r,s,bE] = polyfilt(d,m,x,t,[],[],gcf);  % Takto prikreslime
  title(['SINUSOVKA FILTROVANA S OKNEM DELKY ',num2str(m)]);
  ax=axis;
  hold off;
pause;
%---------------------------------------------------------------
% POKUD POZADUJEME PREDIKCI, STACI ZADAT VHODNE VYSTUPNI SIT tt:
% Zkusime predikci o 5 hodnot do minulosti a soucasne o 5 hodnot
% do budoucnosti:
  tt = -5:n-1+5;
  yy = 2 + sin(2*pi*tt./n);  % Idealni sinusovy trend
  figure; plot(tt,yy,'r','LineWidth',2); hold on
  [ys,E,r,s,bE] = polyfilt(d,m,x,t,tt,[],gcf);  % Takto prikreslime
  title(['PREDIKOVANA SINUSOVKA FILTROVANA S OKNEM DELKY ',num2str(m)]);
% Pouzijeme stejny rozsah os y jako v predchozim grafu
  ax(1:2)=[min(tt),max(tt)]; axis(ax);
  hold off;
pause;
%---------------------------------------------------------------
% PREDCHOZI POSTUP NYNI ZOPAKUJEME AVSAK NA 10x JEMNEJSI VYSTUPNI SITI tt;
% BUDEME TEDY VLASTNE PROVADET I INTERPOLACI MEZI BODY x:
  tt = -5:0.1:n-1+5;
  yy = 2 + sin(2*pi*tt./n);  % Idealni sinusovy trend
  figure; plot(tt,yy,'r','LineWidth',2); hold on
  [ys1,E1,r1,s1,bE1] = polyfilt(d,m,x,t,tt,[],gcf);  % Takto prikreslime
  title(['PREDIKOVANA SINUSOVKA FILTROVANA S OKNEM DELKY ',num2str(m)]);
% Pouzijeme stejny rozsah os y jako v predchozim grafu
  ax(1:2)=[min(tt),max(tt)]; axis(ax);
  hold off;
  zoom on
% NYNI JE ZAPNUT REZIM ZVETSOVANI.
% Kazde kliknuti mysi levym/pravym tlacitkem v miste vyhlazene krivky
% zvetsi/zmensi detail v jeho okoli.
% Vsimnete si nespojitosti vyhlazene krivky zpusobenych tim, ze lokalni
% regresni polynomy na sebe spojite nenavazuji. Tento efekt odstranuje
% pouziti tzv. splajnu, coz jsou funkce spojite slepene z polynomu.
keyboard;
  zoom off
%---------------------------------------------------------------
% UKAZEME NYNI JAK LZE PREDCHOZI VYPOCET EKVIVALENTNE ZOPAKOVAT
% POMOCI OPERATOROVE METODY:
  [FIR,E0] = polyfilt(d,m,[],t,tt);  % FIR=operatorova matice
                                     % E0=[E0s,E0p]=konfid.int. pri sigma=1
% FIR nezavisi na datech, takze ji lze pozdeji opakovane pouzit
% pro vyhlazeni libovolnych dat x definovanych na siti t s vyhlazenym
% vystupem ys2 na siti tt:
  ys2=FIR*x;  % Vyhlazeni
  E2 =s*E0;   % Skutecne konfidencni intervaly dopocteme prenasobenim
              % odhadem sigma (smerodatna odchylka chybove slozky);
              % zde jsme pouzili hodnotu s z predchoziho vypoctu;
  max(abs(ys1-ys2))
  max(abs(E1-E2))
% Vidime, ze vyhlazene prubehy i konfidencni intervaly jsou totozne
% v ramci numericke presnosti.
pause;
%---------------------------------------------------------------
% POSTUPNA POLYNOMIALNI REGRESE PREJDE V OBYCEJNOU REGRESI, POKUD
% delka filtru je nejvetsi mozna, tj. je presne rovna poctu dat: m=n
  figure; plot(tt,yy,'r','LineWidth',2); hold on
  [ys1,E1,r1,s1,bE1] = polyfilt(d,n,x,t,tt,[],gcf);
  title(['PREDIKOVANA SINUSOVKA FILTROVANA S OKNEM DELKY ',num2str(n)]);
%  ax(1:2)=[min(tt),max(tt)]; axis(ax);
  hold off;
pause;
%---------------------------------------------------------------
% Tentyz vysledek bychom meli obdrzet pomoci 'polfs', pokud do modelu
% nezadame zadne harmonicke komponenty (kh=[]):
  [ysE2,trE2,pE2,r2,s2,bE2] = polfs(x,d,[],t,tt,[],figure);
  all(pE2==0)             % Chybi periodicka komponenta =>
  max(abs(ysE2-trE2))     % Vyhlazeny prubeh = trend
% Porovnani:
  max(abs(ys1-ysE2(:,1)))
  max(abs(E1 -ysE2(:,2:3)))
  s1,s2                   % Dobry odhad sigma=0.2 pouzite pri generovani
  max(abs(bE1-bE2))
pause;
%---------------------------------------------------------------
% Pro ekvivalentni vypocet lze take pouzit zabudovane prikazy
% 'polyfit' a 'polyconf' (jen statisticky toolbox):
  t0=(t(1)+t(n))./2;
  [p,S] = polyfit(t-t0,x,d);              % Vypocte regresni polynom
  [ys3,E3] = polyconf(p,tt.'-t0,S,0.05);  % Vypocte JEN predikcni
                                          % konfidencni interval
% Porovnani:
  max(abs(ys1-ys3))
  max(abs(E1(:,2)-E3))
  max(abs(bE1(:,1)-p.'))
%================================== KONEC ==================================
echo off
