%% Signálový processing a procesory
%% Anotace
% Zpracování signálů.
% Systémy reálného času.
% Univerzální a signálový procesor.
% Procesory dSPACE a LabView.
% Tvorba programu pro zpracování signálu, RTW.
% Signál, druhy signálů, popis diskrétního signálu a diskrétního systému, obecný postup zpracování signálu.
% Shannon-Kotělnikovův vzorkovací teorém.
% Fitrace, druhy filtrů a jejich popis.
% Tradiční spojité filtry.
% Návrh nerekurzivních filtrů.
% Návrh rekurzivních filtrů.
% Nestandardní filtrace, zpracování signálů, antialiasing.
% Umělé vidění, pořízení obrazů, zpracování obrazů, reprezentace obrazů.
% Předzpracování obrazů.
% Základy rozpoznávání obrazů.
% Čidla založená na umělém vidění.
%
%% Zpracování signálů
% Zpracování signálů.
% Shannon-Kotělnikovův vzorkovací teorém
% Signál, druhy signálů, popis diskrétního signálu a diskrétního systému, obecný postup zpracování signálu.
%% Real-time
% Systémy reálného času.
% Univerzální a signálový procesor.
% Procesory dSPACE a LabView.
% Tvorba programu pro zpracování signálu, RTW.
%% Filtrace
% Fitrace, druhy filtrů a jejich popis.
% Tradiční spojité filtry.
% Návrh nerekurzivních filtrů.
% Návrh rekurzivních filtrů.
% Nestandardní filtrace, zpracování signálů, antialiasing.
%% Zpracování obrazu
% Umělé vidění, pořízení obrazů, zpracování obrazů, reprezentace obrazů.
% Předzpracování obrazů.
% Základy rozpoznávání obrazů.
% Čidla založená na umělém vidění.
%
%% Filtry
% elektrické obvody měnící spektrum vstupních signálů,
% napěťový přenos: P = \frac{u_{1}}{u_{1}} = \frac{u_{out}}{u_{in}},
% amplitudová charakteristika: A(\omega)_{dB} = 20 \log P
% mezní kmitočet: f_c ... pokles o 3 dB
%% Spojité filtry
% Elektronika.
%% Pasivní analogový
% doc. Papežová -- Elektronika pro aplikovanou mechaniku
% Dolní propust:
%
% přenos P(\omega) = \frac{1}{1+j\omega C R},
% amplitudová charakteristika A(\omega) = \frac{1}{\sqrt{1 + (\omega C R)^2}
% fázová charakteristika \varphi(\omega) = {\rm arctan}(-\omega R C)
%
% Horní propust
%
% přenos P(\omega) = \frac{j\omega C R}{1+j\omega C R},
% amplitudová charakteristika A(\omega) = \frac{\omega C R}{\sqrt{1 + (\omega C R)^2}
% fázová charakteristika \varphi(\omega) = {\rm arctan}(\frac{1}{j\omega C})
%
% Filtry vyšších řádů
% Vyšší strmost,
% složitější obvody
%% Aktivní filtry
% Využití operační zesilovačů:
% usnadňují návrh a praktickou realizaci účinných filtrů,
% kompenzují útlum způsobený pasivními součástkami,
% umožňují dokonalé oddělení vstupu a výstupu filtru,
% umožňují impedanční přizpůsobení,
% lze je řadit do kaskády,
% realizace filtru bez použití indukčnosti
% Dolní propust (lowpass)
%
% Horní propust (hipass)
%
%
%% Návrh filtru
% Nejjednodušší pomocí odvozených vztahů pro normované filtry, přesně platí pouze pro ideální prvky,
% nástroje: MATLAB, Texas Instruments Filter Designer, ...
%% Analog vs Digital
% Analogové
% Levné,
% rychlé,
% velký dynamický rozsah.
% nevhodné pro nízké frekvence,
% nepřesné,
% závislé na vlivech okolí (teplota, ...)
% Digitální
% Přesné,
% univerzální,
% vhodné i pro nízké frekvence,
% stabilní.
%% Digitální filtry
%% FIR (Finite Impulse Response)
% Zaručená stabilita, výstup závislý pouze na vstupech do filtru.
% y_{\rm FIR}(k) = f(u(k), \ldots, u(k-N)) = w_1u(k) + w_2 u(k-1) + \ldots
% Výstup zpožděn o skupinové zpoždění filtruM = \frac{N-1}{2} ,
% lineární fázová charakteristika.
%% IIR (Infinite Impulse Response)
% Nemají zaručenou stabilitu, výstup závislý na vstupech a minulých výstupech.
% y_{\rm IIR}(k) = f(u(k), \ldots, u(k-N_u), y(k-1), \ldots , y(k-N_y)
% Dosahují stejných vlastností jako FIR při nižších řádech, tedy
% nižší výpočetní složitost,
% kratší zpoždění mezi vstupem a výstupem.
% Nelineární fáze, obecný návrh je složitější oproti FIR.
% Citlivé na numerickou přenost.
%% Příklady
% Diskrétní integrátor:
% Obdélníkové pravidlo: y(k) = y(k-1) + x(k)
% Lichoběžníkovo pravidlo: y(k) = y(n-1) + \frac{1}{2} \big(x(k) + x(k-1)\big)
% Diskrétní derivátor:
% y(k) = x(k) - x(k-1)
% y(k) = \frac{1}{2} \big(x(k+1)-x(k-1)\big)
%% Praktická aplikace
%% Offline použití
% Použiju nejlepší, jaký mám -- možno použít nekauzální filtr.
%% Online použití
% 🤷🏻‍♂️🤷🏻‍♂️
% Podle aplikace.

filterDesigner

%
%% Adaptivní filtr
% Koeficienty filtru se mění na základě optimalizačního kritéria za běhu procesu.
%% Aplikace
% Použití HONU/LNU (high order neural unit / linear neural unit) k identifikaci kmitavého systému
% Online identification using linear neural unit with guaranteed weights convergence
% Systém

plant = tf(50000*[0.6106    0.0753e3],[1       204.21        15278   4.1545e+05   5.8866e+05]);

%
% "Měření"

fS = 100;
t = 0:1/fS:20;
u = .5*square(t*2*pi); % 2*pi*t = 1 Hz
u(fix(numel(u)/2):end) = .4*u(fix(numel(u)/2):end)+0.3*square(t(fix(numel(u)/2):end)*2*pi*3.468135); 
y = lsim(plant,u,t);

figure
tiledlayout(2,1,"Padding","compact","TileSpacing","compact")
nexttile
plot(t,u,"LineWidth",1.3)
ylabel("u")
grid on
nexttile
plot(t,[y],"LineWidth",1.3)
ylabel('y')
grid on
xlabel('t [s]')

%
% Parametry neuronu

nU = 10;
nY = 2;
LNUmaxN = nU+nY;

mu = 2;
epsilon = 1e-5;
wLNUinit = 0*randn(1,1+2*LNUmaxN);
wBuff = zeros(2*LNUmaxN+1);
wChange = zeros(size(y));
yLNUhist = zeros(size(y));
eHist = zeros(size(y));
SSEHist = zeros(size(y));
rhoHist = zeros(size(y));
muHist = zeros(size(y));
eHistTemp = zeros(LNUmaxN,1);

%
% Simulace online běhu

wPrev = wLNUinit;
k0 = max([find(t>1,1),LNUmaxN+1]);
muHist(1:k0)=mu;
SSEHist(k0) = sum(y(k0-LNUmaxN:k0).^2);
for Kdx = 2:numel(t)
    if Kdx<LNUmaxN+1
        uVec= zeros(LNUmaxN,1);
        uVec(1:Kdx) = u(Kdx:-1:1);
        yVec = zeros(LNUmaxN,1);
        yVec(1:Kdx) = y(Kdx:-1:1);
    else
        uVec = u(Kdx:-1:Kdx-LNUmaxN+1); % u up to k
        yVec = y(Kdx:-1:Kdx-LNUmaxN+1); % y up to k-1, (otherwise might lead to yLNU(k) = y(k))
    end
    
    [yLNUhist(Kdx),wPrev,rhoHist(Kdx),muHist(Kdx)] = LNUfcn(1,uVec,yVec,nU,nY,muHist(Kdx-1),epsilon,wPrev,SSEHist(Kdx-1),LNUmaxN);
    wBuff = circshift(wBuff,1,1);
    wBuff(1,:) = wPrev;
    wChange(Kdx) = mean(abs(wBuff(2,:)-wBuff(1,:)));
    eHist(Kdx) = y(Kdx) - yLNUhist(Kdx);
    if Kdx>k0
    
    end
    if Kdx<1+LNUmaxN
        eHistTemp(1:Kdx) = eHist(1:Kdx);
        eHistTemp(Kdx+1:end) = max(abs(eHist));
        SSEHist(Kdx) = sum(eHistTemp.^2);
        
    else
        SSEHist(Kdx) = sum(eHist(Kdx-LNUmaxN:Kdx).^2);
        
    end
end
SSEHist(1:10)=SSEHist(11);

figure
tiledlayout(4,1,"Padding","compact","TileSpacing","compact")
nexttile([3,1])
plot(t,[yLNUhist],"LineWidth",1.3)
hold on
plot(t,[y],'--',"LineWidth",1.3)
legend({"y_m","y"})
grid on
title(sprintf("n_u = %i, n_y = %i, mu = %.3g, SSE_{min} = %.2g",nU,nY,muHist(end),min(SSEHist)))
nexttile
semilogy(t,SSEHist,"LineWidth",1.3)
grid on
ylabel("SSE")

%
% Funkce zajišťující adaptaci neuronu

function  [yLNU,w,rho,mu] = LNUfcn(learn,uVec,yVec,nU,nY,mu,epsilon,wPrev,SSE,maxN)
LNUmaxN = maxN;

colX = zeros(1,1+2*maxN);
colX(1)=1;
yLNU = 0;
yLNU = yLNU+1*wPrev(1);
for idx=1:nU
    yLNU=yLNU+wPrev(idx+1)*uVec(idx);
    colX(idx+1)=uVec(idx);
end
for idx=1:nY
    yLNU=yLNU+wPrev(idx+1+LNUmaxN)*yVec(idx+1); %y(k) is only for e(k)
    colX(idx+1+LNUmaxN)=yVec(idx);
end

mun = mu/(epsilon+norm(colX.*colX));
rho = StabEval(colX([1 2:nU+1 LNUmaxN+1+1:LNUmaxN+1+nY]),mun);
if round(rho,8)>1
    mu=0.9*mu;
    w=wPrev;
    return
else
    if rho<1
        mu=1.01*mu;
    end
end

% if learning is disabled
if learn<1 
    w=wPrev;
    return
else
    % weights update
    e = yVec(1) - yLNU;
    w = wPrev+mun*e*colX;
end

end

function rho = StabEval(x,mun)
    Amat = eye(numel(x))-mun*(x'*x);
    rho = norm(Amat,2); % frobenius norm
    rhoS = max(abs(eig(Amat))); % spectral radius
    if rho~=rhoS
        sprintf("spectral radius: %f, L2 norm: %f",rhoS,rho);
    end
end