% Nada Amin (namin@mit.edu)
% 6.003 Lab 1
% Matlab Code
load lineup;
sound(y, 8192);
N = 1000;
alpha = 0.5;
% (a)
a = [1 zeros(1, N)];
b = [1 zeros(1, N-1) alpha];
d = [1 zeros(1, N)];
he = filter(b, a, d);
figure;
plot([0:N], he);
xlabel('n'); ylabel('he[n]');
% (b) analytical only
% (c)
her = filter(a, b, [1 zeros(1,4000)]);
figure;
plot([0:4000], her);
xlabel('n'); ylabel('her[n]');
% (d)
z = filter(a, b, y);
figure;
plot([0:length(z)-1],z);
xlabel('n'); ylabel('z[n]');
sound(z, 8192);
% (e)
hoa = conv(he, her);
figure;
plot([0:length(hoa)-1],hoa);
xlabel('n'); ylabel('hoa[n]');
% (f)
[Ny1, Ryy] = detect1(y);
half = (length(Ryy)-1)/2;
figure;
plot([-half:half], Ryy);
xlabel('n'); ylabel('Ryy[n]');
ay1 = alpha1(z, y, Ny1);

[Ny2, Ry2y2] = detect1(y2);
half = (length(Ry2y2)-1)/2;
figure;
plot([-half:half], Ry2y2);
xlabel('n'); ylabel('Ry2y2[n]');
ay2 = alpha1(z, y2, Ny2);

% y3 is composed of two echoes
% will process it in two passes
% **pass 1**
[N1, Ry3y3] = detect1(y3);
half = (length(Ry3y3)-1)/2;
figure;
plot([-half:half], Ry3y3);
xlabel('n'); ylabel('Ry3y3[n]');
% approximate alpha to be the ratio of the current little peak with the
% original center peak
ci = half+1;
Rxx = conv(z, flipud(z));
a1 = Ry3y3(ci-N1)/Rxx(ci);
% round alpha to the nearest decimal
a1 = round(a1*10)/10;
% z1 is y3 filtered for this echo
z1 = filter([1 zeros(1, N1)], [1 zeros(1, N1-1) a1], y3);
% **pass 2**
[N2, Rz1z1] = detect1(z1);
half = (length(Rz1z1)-1)/2;
figure;
plot([-half:half], Rz1z1);
xlabel('n'); ylabel('Rz1z1[n]');
L = length(z);
if mod(L,2)==0
    c = L/2;
else
    c = (L+1)/2;
end
n = c+max(N2,N1)+1;
a2 = (y3(n)-z(n)-a1*z(n-N1))/z(n-N2);

% sanity checks
zf2 = filter([1 zeros(1, Ny2)], [1 zeros(1, Ny2-1) ay2], y2);
dzf2 = sum(abs(z-zf2)); % must be very small

if N1 < N2
    minN = N1;
    maxN = N2;
    mina = a1;
    maxa = a2;
else
    minN = N2;
    maxN = N1;
    mina = a2;
    maxa = a1;
end

as = [1 zeros(1, maxN)];
bs = [1 zeros(1, minN-1) mina zeros(1, maxN-minN-1) maxa];
zf3 = filter(as, bs, y3);
dzf3 = sum(abs(z-zf3)); % must be very small
