% 6.003 Problem Set 3, Problem 7
% WonRon Cho (leopardi)
% Nada Amin (namin)

N = 8;
y = [0 -1 -1 -1 1 1 0 1];
k = [0:N-1];
n = [0:N-1];
% find the fourier series of y[n] (b_0 to b_7)
b = zeros(1,N);
for k=0:N-1,
  b(k+1) = (1/N)*sum(y .* (cos(-2*pi*k*n/N) + j*sin(-2*pi*k*n/N)));
end
% b =  [0            -0.1250 + 0.4786i   0.2500 + 0.0000i  -0.1250 + 0.2286i        0 - 0.0000i  -0.1250 - 0.2286i   0.2500 + 0.0000i  -0.1250 - 0.4786i]
       
% now let's find a_k by using part y[n] = x[n] - x[n-1]
% a_k = b_k / (1 - (cos(-2*pi*k/N) + j*sin(-2*pi*k/N)))
k = [0:N-1];
a = b ./ (1 - (cos(-2*pi*k/N) + j*sin(-2*pi*k/N)));
% a [NaN             0.5152 + 0.3902i   0.1250 - 0.1250i  -0.0152 + 0.1402i  -0.0000 - 0.0000i  -0.0152 - 0.1402i   0.1250 + 0.1250i   0.5152 - 0.3902i]


% the problem is a_0
% we're going to use fact b) and Parseval's Relation
abs_ao = sqrt(12/N - sum(abs(a([2:N])).^(2*ones(1,N-1))));
% abs_ao = 0.75


% we want neat cos and sin for our answer
%cospart = acos(real(a))/pi
%sinpart = acos(imag(a))/pi

% now we use the synthesis equation
almost_x_0 = sum(a([2:N]));
%almost_x_0 =  1.2500 + 0.0000i

% ok, so a_0 could be either 0.75 or -0.75
% but we can use the difference equation
% y[1] = x[1] - x[0]
% x[1] = x[0] + y[1]
% if a_0=0.75, we have x[1] = 1.25+0.75-1 = 1
% if a_0=-0.75, we have x[1] = 0.5-1=-0.5
% sigh, they both work
a(1) = 0.75;

% finally, a_k
% a = [0.75 0.5152+0.3902i   0.1250-0.1250i  -0.0152+0.1402i  0  -0.0152-0.1402i  0.1250+0.1250i   0.5152-0.3902i]

% checks
yc = [0:N-1];
for k=0:N-1,
  yc(k+1) = sum(b .* (cos(k*2*pi*n/N) + j*sin(k*2*pi*n/N)));
end

xp = [0:N-1];
for k=0:N-1,
  xp(k+1) = sum(a .* (cos(k*2*pi*n/N) + j*sin(k*2*pi*n/N)));
end

a(1) = -0.75;
xn = [0:N-1];
for k=0:N-1,
  xn(k+1) = sum(a .* (cos(k*2*pi*n/N) + j*sin(k*2*pi*n/N)));
end

en = sum(abs(xn).*abs(xn))
ep = sum(abs(xp).*abs(xp))

% they both give 12!
% so it seems like both a_0 is 0.75 and a_0 is -0.75 are ok???

yp = zeros(1,N);
yp([2:N]) = xp([2:N]) - xp([1:N-1]);
yp(1) = xp(1) - xp(N);

yn = zeros(1,N);
yn([2:N]) = xn([2:N]) - xn([1:N-1]);
yn(1) = xn(1) - xn(N);


sum(abs(y-yp))
sum(abs(y-yn)) % both are very small
