% ECE4703B10
% This code allows you to do system identification using the brute-force
% correlation + matrix inversion method
% DRB 09-Dec-2010

[s,fs] = wavread('lab6part1_sys1.wav'); % load the input/output of the unknown system
N = 20;                                 % set filter length

x = s(:,1);                             % input to unknown system
d = s(:,2);                             % output of unknown system

% compute the p vector
tmp = xcorr(d,x,N-1,'unbiased');
p = tmp(N:end);

% compute the R matrix
tmp = xcorr(x,x,N-1,'unbiased');    % compute r_k,ls
r = tmp(N:end);
R = zeros(N,N);
for k=1:N,
    for ell=1:N,
        R(k,ell) = r(abs(k-ell)+1);
    end
end

bmmse = inv(R)*p;
stem(0:N-1,bmmse);
xlabel('filter coefficient');
ylabel('value of filter coefficient');
grid on