LMS(FIR_A,FIR_B,L,wave)

function [w_out mse_out ref_out] = LMS(FIR_A,FIR_B,L,wave)
% [w mse ref res iter] = LMS(FIR_A,FIR_B,L,wave)
% LMS filter to solve the system identification problem represented below:
% 
%            ----------               ----------
%     ------- FILTER A -----out_A----- FILTER X ---out--
%     |      ----------               ----------       |
% in  |                                                |
% ----|                                                |+
%     |      ----------                            - -----
%     ------- FILTER B -----out_B-------------------- SUM ---error---
%            ----------                              -----
%
% FILTER_X is unknown and to be derived. This problem is called "filter matching"
% and is encountered when one needs to augment a certain filter (A)
% in order to match the behavior of a reference filter (B).
%
% Parameters have the following significance:
%   FIR_A:  impulse response (taps) of FILTER_A as a column vector
%   FIR_B:  impulse response of FILTER_B as a column vector
%   L:      number of weights in the unknown filter
%   wave:   reference waveform. Available options: 
%           ('square', 'triangle', 'noise', 'chirp', 'multiburst')
%
% The outputs have the following significance:
%   w_out:      weights of the unknown filter
%   mse_out:    evolution of the mse with iterations
%   ref_out:    input (reference) signal
%-------------------------------------------------------------------------

% constants
mu = 10^-5;         	% learning rate

% parameters for the reference signal
n_samples = 20000;  	% number of samples in the reference signal
periods = 10;        	% number of periods in the input square signal
signal_amplitude = 10;	% amplitude of the reference signal
lead_border = 500;   	% skipped samples at the beginning
lag_border = 500;     	% skipped samples in the end
samples_per_period = int32(n_samples / periods);

% parameters specific to the square wave
duty_cycle = 0.7;     	% duty-cycle for the square wave input

% parameters specific to the triangular wave
skewness = 0.7;       	% skewness of the triangular wave

% parameters specific to the periodic chirp signal
T0 = samples_per_period / 4;   % initial period in samples
Tf = samples_per_period / 40;  % final period in samples

% parameters specific to the multiburst signal
flat_line = samples_per_period / 5; % samples spent on flat region that separates the bursts

% preliminary steps
fir_A = double(FIR_A);
fir_B = double(FIR_B);
[len_Ay len_Ax] = size(fir_A);
[len_By len_Bx] = size(fir_B);  
len_A = max(len_Ax,len_Ay);
len_B = max(len_Bx,len_By);
if len_Ax > len_Ay
    fir_A = fir_A';
end
if len_Bx > len_By
    fir_B = fir_B';
end
max_len = max(len_A, len_B);

% normalize filters
fir_A = double(fir_A)/double(sum(fir_A));
fir_B = double(fir_B)/double(sum(fir_B));

% generate the input reference signal. 
in = zeros(n_samples,1);

% avoid division by zero in the triangular wave case
if skewness == 1 
    skewness = 1.0-1.0/samples_per_period;
end

if strcmp(wave,'square')
    samples_per_high = int32(samples_per_period * duty_cycle);
    for i=1:n_samples
        local_idx = mod(i,samples_per_period);
        if local_idx < samples_per_high
            in(i)= signal_amplitude;
        end
    end
elseif strcmp(wave,'triangle')
    samples_per_ascending_front = samples_per_period * skewness;
    for i=1:n_samples
        local_idx = mod(i,samples_per_period);
        if local_idx < samples_per_ascending_front
            in(i)= signal_amplitude * double(local_idx) / ... 
                double(samples_per_ascending_front);
        else 
            in(i) = double(signal_amplitude) * double(samples_per_period - local_idx)/ ... 
                double(samples_per_period - samples_per_ascending_front);
        end
    end
elseif strcmp(wave,'noise')
    in = signal_amplitude * randn(n_samples,1);
elseif strcmp (wave, 'chirp')
    for i=1:n_samples
        local_idx = mod(i,samples_per_period);
        T = T0 - double(local_idx) * double(T0 - Tf)/double(samples_per_period);
        in(i) = double(signal_amplitude) * sin(2 * pi / double(T) * double(local_idx));  
    end
else
    display('Unknown waveform option');
end

% generate signals; This step can be replaced with canned convolution
% but one must be careful about misalignment when filters have unequal lengths
out_A = zeros(n_samples,1);
out_B = zeros(n_samples,1);

for i=1:n_samples
    s = 0;
	for j=1:len_A
        signal_idx = min(n_samples,max(1,i+j-1-int32(len_A/2)));
        kernel_idx = j;
        s = s + in(signal_idx) * fir_A(kernel_idx);
    end
    out_A(i,1)=s;
end

for i=1:n_samples
    s = 0;
	for j=1:len_B
        signal_idx = min(n_samples,max(1,i+j-1-int32(len_B/2)));
        kernel_idx = j;
        s = s + in(signal_idx) * fir_B(kernel_idx);
    end
    out_B(i,1)=s;
end

% make sure we leave enough samples near the borders
lead_border = max( lead_border, L);
lag_border = max(lag_border, L);

% initialize iterations
w = zeros(L,1);
mse = zeros(n_samples,1);
iter = 0;
% main iteration of the LMS algo
for k=lead_border:n_samples - lag_border
    iter = iter + 1;
    y = 0;
    for i=1:L
        y = y + out_A(k+i) * w(i);
    end
    error = out_B(k+int32(L/2+1)) - y;
    mse(k) = error * error;
    lr = 2 * mu * error;
    
    % weight update
    for i=1:L
        w(i) = w(i) + lr * out_A(k+i);
    end 
end

w_out = w;
ref_out = in;
mse_out = mse;