% -----------------------------------------------------------
% System identification via LMS example code
% ECE4703
% DRB 08-Dec-2011
% -----------------------------------------------------------
% User parameters
% -----------------------------------------------------------
fs = 44100;             % sampling frequency
L = 60;                 % length of our adaptive filter
mu = 1/L;               % step size
T = 0.25;               % duration of simulation
iterations = 100;         % iterations for averaging squared error
% -----------------------------------------------------------

% PRELIMINARY CALCULATIONS
Ts = 1/fs;              % sampling period
t = 0:Ts:T;             % time vector
N = length(t);          % length of time vector
hhat = zeros(L,N);      % initialize adaptive filter
y = zeros(1,N);         % output of adaptive filter
e = zeros(1,N);         % error

% CREATE "UNKNOWN" SYSTEM
Fpass = 3000;            % Passband Frequency
Fstop = 4000;            % Stopband Frequency
Dpass = 0.057501127785;  % Passband Ripple
Dstop = 0.031622776602;  % Stopband Attenuation
dens  = 20;              % Density Factor
[M, Fo, Ao, W] = firpmord([Fpass, Fstop]/(fs/2), [1 0], [Dpass, Dstop]);
h  = firpm(M, Fo, Ao, W, {dens});
a = 1;

mse = zeros(1,N);

tic

for i=1:iterations,
    
    % CREATE WHITE NOISE INPUT SIGNAL ON [-1,1]
    x = 2*rand(1,N)-1;
    
    % COMPUTE OUTPUT OF "UNKNOWN" SYSTEM
    d = filter(h,a,x);
    
    % MAIN LMS LOOP
    for n=L:N-1,
        xx = x(n:-1:n-L+1)';        % past input buffer (vectorized)
        y(n) = hhat(:,n)'*xx;       % output of adaptive filter (vectorized)
        e(n) = d(n) - y(n);         % error
        hhat(:,n+1) = hhat(:,n) + mu*e(n)*xx;   % update coefficients (vectorized)
    end
    
    % ACCUMULATE STATISTICS FOR MSE PLOT
    mse = mse+e.^2/iterations;
    
    [i toc]

end

% MAKE PLOTS
figure(1)
subplot(3,1,1)
stem(0:L-1,hhat(:,L));
hold on
stem(0:M,h,'r');
hold off
legend('adaptive filter','unknown system');
title('initial');
subplot(3,1,2)
stem(0:L-1,hhat(:,L+200));
hold on
stem(0:M,h,'r');
hold off
legend('adaptive filter','unknown system');
title('after 200 samples');
subplot(3,1,3)
stem(0:L-1,hhat(:,end));
hold on
stem(0:M,h,'r');
hold off
legend('adaptive filter','unknown system');
title('at end of simulation');
xlabel('filter coefficient');

figure(2)
semilogy(t,mse);
xlabel('time');
ylabel('mean squared error');
grid on
