% E05_1  Adaptive Filters: System_Ident      gr_08may20  25apr22
%
% With LP filtered input noise: reduce mu; smaller band --> larger error
%   in the model (higher input correl.). Try e.g. alpha=.5, .8, .9, .99
% If x = sine (deterministic, single freq.), LMS finds a solution for 
%   _that_ input, error goes down but model is totally wrong
% Try also noisy plant
% Try also IIR plant, mmse FIR approx is obtained


clear all;

plant = [1 -.3 .4 -.5];  % unknown system we want to identify
Lp = length(plant);

L = 7;  % length of model (FIR adaptive filter)

N = L * 111;  % N = multiple of L to simplify calc. of Rx
n = 1:N;
 x = randn(1,N); % input signal: white Gaussian noise
% x = filter(1,[1 -.5],x); %  correlated noise. Try e.g. alpha=.5, .8, .9
% x = sin(pi/23*n);  % sine function: NOTE error goes down but model is wrong!

x = x ./ sqrt(mean(x.^2));  % global energy of signal set to 1

Rx = corrcoef(reshape(x,111,N/111))  % input corr. matrix
lambda = sort(eig(Rx),'descend')'   % and its eigenvalues

mu = 0.1;  % mu < 2/lambda_max
% (e.g. mu=0.2 for white noise, mu=0.1 for colored noise, 0.05 for noisy plant)
h_est = [1, zeros(1,L-1)];   %  model init.
for k=L:N
    d(k) = x(k:-1:k-Lp+1) * plant';
 %   d(k) = x(k:-1:k-Lp+1) * plant' + (rand(1,1)-.5)/5;    % noisy plant
 %   a=.7; d(k) = a*x(k) - (1-a)*d(k-1);    %  a different, very simple IIR plant
    y(k) = x(k:-1:k-L+1) * h_est';
    e(k) = d(k) - y(k);
    h_est = h_est + mu * e(k) * x(k:-1:k-L+1);    
end
plant, h_est

qer = e.^2;
subplot(211); plot(n,d,n,y)
subplot(212); semilogy(n,qer)
MSE_Last100 = mean(qer(N-99:N))
MSE_model = mean(h_est - [plant, zeros(1,L-Lp)]).^2

% plant_FIRapprox = impz(a,[1  1-a],L)'  % for case of simple IIR plant

return
