Norm MPEG-1 Layer III (Mp3) Filter banks distortion
I have tested your code. 

There are problems with the decimation and upsample stages. When you run the analysis and synthesis filters at the undecimated rate, then you get the perfect reconstruction (plus a filter shift)

First of all in case of applications which are sensitive to group delay, you should better use the conv / filter function and then compansate for the group delay.

Also, the decimate function applies a lowpass anti-aliasing filter, which kills the signal content at the $k$-th channel considering the fact that its output is a bandpass signal rather than lowpass. Hence you should avoid the anti-aliased decimation but instead employ a bandpass decimation strategy.
A direct subsampling by 32 would provide a convincing result.

Same goes with interpolation part. After you expand your signal, you should select the proper band of image spectrum by a suitable bandpass filter, rather than a lowpass filter from a conventional interpolator. A direct expansion may convince as its output will already be bandpass filtered later by the polyphase filterbank.

That being said, your main problem is the document that you based your code. That's far too out of signal processing context. I suggest you read the MPEG audio standard book of BOSI. But it's up to you of course.

The following is a working corrected MATLAB / OCTAVEcode:

% S0 - Load the prototype lowpass filter impulse response:
% --------------------------------------------------------
load h2.mat;         % h[n] is the prototype lowpass filter... of length 512
L = length(h);

figure,subplot(2,1,1)
stem([0:L-1],h);title('The Prototype Lowpass Filter');
subplot(2,1,2)
plot(linspace(-1,1,4*L),20*log10(abs(fftshift(fft(h,4*L)))));
grid on;


% S1 - Create the 32 band filter bank by cosine modulation
% --------------------------------------------------------
numbands = 32;                  % number of bands
n=0:L-1;

hha=zeros(numbands,L);         % band of filters
for k=0:1:numbands-1 
   hha(k+1,:) = h.*cos( ( (2*k+1)*pi*(n-16) ) / (2*numbands) );
end

figure
plot(linspace(-1,1,4*L),20*log10(abs(fftshift(fft(hha(1,:),4*L)))));
title('ANALYSIS FILTER BANK');
hold on
for k=1:numbands
    plot(linspace(-1,1,4*L),20*log10(abs(fftshift(fft(hha(k,:),4*L)))));
end

% S2 - Design the synthesis filterbank:
% -------------------------------------
numbands = 32;                     % number of banks
n=0:L-1;

hhs = zeros(numbands,L);           % band of filters
for k=0:1:numbands-1 
   hhs(k+1,:) = h.*cos( ( (2*k+1)*pi*(n+16) ) / (2*numbands) );
end


% S3 - Generate a test input signal:
% ----------------------------------
N = 2*1024;      % NOTE: According to the MPEG specification a determined LENGTH must be applied.
wav_in = sin(0.05*pi*[0:N-1]);        % pure sine tone

figure,subplot(2,1,1)
plot(wav_in);title('input signal')
subplot(2,1,2)
plot(linspace(-1,1,4*N),20*log10(abs(fftshift(fft(wav_in,4*N)))));


% S4 - Apply a wav file to it,  get every output of the every sub-band:
% --------------------------------------------------------------------
yyd = zeros( numbands, floor(N/numbands));   % decimated outputs..
yyu = zeros( numbands , N);

for k=1:1:numbands
    temp = conv( wav_in,hha(k,:));
    yyu(k,:) = temp(L/2+1:L/2+N);
    yyd(k,:) = yyu(k,1:numbands:end);
end


% S6 - Apply synthesis on the decimated signal
% --------------------------------------------
ys = zeros(1, N);
yi = zeros(1,N);

for k=1:1:numbands
    yi(1:numbands:end) = yyd(k,:);
    temp = conv(yi,hhs(k,:));
    ys = ys + temp(L/2+1:L/2+N);
end
ys = numbands*ys;

% Display the result:    
figure,subplot(2,1,1)
plot(ys);title('Synthesized Back');
subplot(2,1,2)
plot(linspace(-1,1,4*N),20*log10(abs(fftshift(fft(ys,4*N)))));


With the outputs being:

1- The prototype filter:



2- 32 channel filter bank:



3- A pure tone input:



4- An reconstructed back output: