%% Matérn and squared-exponential kernels: intuition from a control/signal processing engineer point of view
% © 2024, Antonio Sala, Universitat Politècnica de València, Spain. All rights reserved.
% Video presentation at: https://personales.upv.es/asala/YT/V/maternEN.html
% This code was run with Matlab R2024a
% Objectives: understand the "filtering" intuition (spectral factor transfer function) on Matérn kernel and the squared-exponential one.
%
%% Frequency response
% The "Matérn" kernel [ from Bertil Matérn (1917 – 2007), Wikipedia ] is a generalisation to maybe "fractional n" of white noise filtered by \frac{A}{(\tau s+1)^n} , well actually more similar to this \frac{A}{(\frac{\rho}{\sqrt{2n-1}} s+1)^n}... and in fact, A also changes with n so that the 2-norm (steady-state standard deviation) keeps constant to a pre-defined value \sigma.
% But these issues involving Bessel functions, Gamma functions, etc. are not needed in detail in order to intuitively understand what we are doing if we come from a control systems background...
%
% So, apart from "amplitude correction" to keep variance constant, we are doing something like this:

s=tf('s');
ntst=[1 2 3 6 30];
Gtst=cell(5,1);Gtst2=cell(5,1);
for k=1:5
    n=ntst(k);
    tau=1/sqrt(2*n-1);
    Gtst{k}=1/(tau*s+1)^n;
    Gtst2{k}=Gtst{k}/norm(Gtst{k});
end
w_tst=logspace(-1.5,1);
bodemag(Gtst{:},w_tst), grid on, ylim([-20 2])
title("Bode diagram (Spectral Factor),   PSD -> multiply dB*2")
hold on
plot(w_tst,20*log10(exp(-w_tst.^2/2/2)),"-.",LineWidth=2),
hold off
legend("1", "2", "3", "6", "30", "\infty")

% This is the actual amplitude correction so the integral of the "power spectral density" (i.e., variance) keeps constant to one:

for k=1:length(Gtst)
    Gtst2{k}=Gtst{k}/norm(Gtst{k});
end
bodemag(Gtst2{:},logspace(-1.5,1)), grid on, ylim([-20 8])

% Let us check Wikipedia formulae...

A=3;tau=1; n=5;

% Causal spectral factor:

G=@(s) A/(tau*s+1)^n;

% Power spectral density:

syms w real
psd=simplify(G(1j*w)'*G(1j*w),60)

% Auto-covariance:

syms x real %distance between points
kk=simplify(ifourier(psd))

% Integral of PSD is the variance, its. square root is standard deviation of the stationary process.

sg=simplify(sqrt(int(psd,-inf,inf))/sqrt(2*sym(pi))); %standard deviation in steady state
eval(sg)
sg1=norm(G(s)) %2-norm of a dynamic system
eval(sqrt(subs(kk,0))) %autocovariance at zero distance is variancell=5+n;

% Now, Matérn kernel formulae from Wikipedia end up in the covariance:

nu=sym(n-1/2);
rho=tau*sqrt(2*nu);
kk2=simplify( sg^2*2^(1-nu)/gamma(nu)*(sqrt(2*nu)*x/rho)^nu*besselk(nu,sqrt(2*nu)*x/rho) )

% *Note that "ifourier" command above produced a polynomial with coefficients of odd degree with |x|, |x|^3 , .... The expression kk2 is supposed to be evaluated with x\geq 0 to yield identical results (it must be used with argument equal to the norm of the difference between two abscissa of the Gaussian process).
% When n\to\infty we get the "squared-exponential" kernel:

expq=exp(-(x/rho)^2/2)*sg^2
ll=5+n; %cosmetic, axis stuff
fplot(kk,[0, ll*tau],LineWidth=4), hold on, grid on
fplot([kk2 expq],[0, ll*tau], LineWidth=2), hold off
title("Autocorrelation")
legend("ifourier from transfer function","bessel,gamma","exp-quadratic")
xlabel("time"),ylabel("covariance")