% Part 1 Codec (wraps around at edges, uses basic averaging for overlap)
%input: vector x of input signal
%       length n of window
%       b quantization bits
%output: vector out of output signal
%Example usage: out=prob1codec(cos((1:2^(12))*2*pi*256/2^(13)), 32, 3);
%               example signal is 1/2 sec. pure tone of frequency f = 256
function out=prob1codec(x, n, b)
len=length(x);
nw=floor(len/n);                                    % number of length n windows in x
x=x(1:n*nw);                                        % cut x to integer number of n windows
Fs=2^(13);                                          % Fs=sampling rate
L=0.3;                                              % [-L,L] amplitude range
q=2*L/(2^b-1);                                      % q used for b bits on interval [-L, L]
for i=1:n                                           % form the MDCT matrix
    for j=1:2*n
        M(i,j)= cos((i-1+1/2)*(j-1+1/2+n/2)*pi/n);
    end
end
M=sqrt(2/n)*M;
N=M';                                               % inverse MDCT

%%%%%%%%%%
% ENCODE %
%%%%%%%%%%
x=0.3*x/max(abs(x));                                % normalize signal to max amplitude = 0.3
out=[];
for k=1:nw-1                                        % loop over length 2n windows
    x0=x(1+(k-1)*n:2*n+(k-1)*n)';
    y0=M*x0;
    y1(:,k)=round(y0/q);                            % transform components quantized
end;
% wrap around last n to first n and append to end of transformation
x_last = cat(2, x(1+(nw-1)*n:nw*n), x(1:n))';
y_last = M*x_last;
y1(:,nw) =round(y_last/q);

%%%%%%%%%%
% DECODE %
%%%%%%%%%%
for k=1:nw-1
    y2=y1(:,k)*q;                                   % transform components dequantizd
    w(:,k)=N*y2;                                    % invert the MDCT
    if k>1
        w2 = w(n+1:2*n,k-1);
        w3 = w(1:n,k);
        %if k>4                                     % plot window for debuggging
        %    plot(1+(k-4)*n:2*n+(k-4)*n, w(:,k-3), ...
        %    1+(k-3)*n:2*n+(k-3)*n, w(:,k-2), ...
        %    1+(k-2)*n:2*n+(k-2)*n, w(:,k-1), ...
        %    1+(k-1)*n:2*n+(k-1)*n, w(:,k  ), ...         
        %    1+(k-4)*n:2*n+(k-1)*n, x(1+(k-4)*n:2*n+(k-1)*n));
        %    pause(5);
        %end
        out=[out;(w2+w3)/2];                        % collect the reconstructed signal
    else
        out=[0*w(1:n,k)];                           % fill first n slots with zeros (to be reconstructed later)
	end
end
out=[out; 0*(1:n)'];                                % fill last n slots with zeros (to be reconstructed later)
% reconstruct "wrapped" window (n from end | n from beginning)
y2=y1(:,nw)*q;
w(:,nw)=N*y2;
w2=w(n+1:2*n,nw);w3=w(1:n,1);                       % calculate first n data points
out(1:n) = (w2+w3)/2;
w2=w(n+1:2*n,nw-1);w3=w(1:n,nw);                    % calculate last n data points
out(n*(nw-1)+1:nw*n) = (w2+w3)/2;

% return decoded signal as horizontal array
out = out';
%threetones = horzcat(x, zeros(1,Fs/4), out, zeros(1,Fs/4), (x-out)*5);
%sound(threetones, Fs);
%pause(length(threetones)/Fs);