clear;
close all;
clc;

N = 1e6;
u_k = randn(1,N);
%u_k = (2*rand(N,1)-1)*sqrt(3);

%% Linear Block

tau = [0,1,2,3,4,5];
y_k1 = zeros(1,N);
y_k1(1) = u_k(1);
for k = 2:N
    y_k1(k) = 0.5*y_k1(k-1) + u_k(k);
end

gamma_dash1 = zeros(1, length(tau));

for i = 1:length(tau)
    sum = 0;
    for k = 1:(N-tau(i))
        sum = sum + u_k(k) * y_k1(k+ tau(i));
    end
    gamma_dash1(i) = sum / (N - tau(i));
end

%% Hammerstein System

miu = zeros(1,N);

for i = 1:N
    miu(i) = 3*u_k(i)+2*u_k(i)^2;
end

y_k2 = zeros(1,N);
y_k2(1) = miu(1);
for k = 2:N
    y_k2(k) = 0.5*y_k2(k-1) + miu(k);
end

gamma_dash2 = zeros(1, length(tau));

for i = 1:length(tau)
    sum = 0;
    for k = 1:(N-tau(i))
        sum = sum + u_k(k) * y_k2(k+ tau(i));
    end
    gamma_dash2(i) = sum / (N - tau(i));
end

%Scale
gamma_dash2 = gamma_dash2/gamma_dash2(1);

%% Wiener System

x_k = zeros(1,N);
x_k(1) = u_k(1);
for k = 2:N
    x_k(k) = 0.5*x_k(k-1) + u_k(k);
end

y_k3 = zeros(1,N);

for i = 1:N
    y_k3(i) = 3*x_k(i)+2*x_k(i)^2;
end

gamma_dash3 = zeros(1, length(tau));

for i = 1:length(tau)
    sum = 0; 
    for k = 1:(N-tau(i))
        sum = sum + u_k(k) * y_k3(k+ tau(i));
    end
    gamma_dash3(i) = sum / (N - tau(i));
end

%Scale
gamma_dash3 = gamma_dash3/gamma_dash3(1);