% Polyphase Quadrature Filterbank: Synthesis
%
% usage:
%   x = pqf_synth(S,h);
% where:
%   S     : matrix of decimated subband outputs (N rows <=> N subbands)
%   h     : coefficients of lowpass prototype filter (cutoff at pi/2N radians)
%   x     : vector of time-domain reconstructed samples
%   
% Copyright 1999 --- Phil Schniter

function u = pqf_synth(S,h);

 N = size(S,1);
 L = size(S,2);
 M = length(h);

 % create synthesis window
 hw = 2*h;		% last sample is zero
 for k=1:2:M/2/N,
   hw(k*2*N+[1:2*N]) = -hw(k*2*N+[1:2*N]);
 end;

 % fast cosine matrix transformation
 Vbar = sqrt(N/2)*dct(S); Vbar(1,:) = sqrt(2)*Vbar(1,:);
 Vnew = zeros(2*N,L);
 Vnew(1:N/2,:) = Vbar(N/2+1:2*N/2,:);
 Vnew(N/2+2:3*N/2,:) = -Vbar(2*N/2:-1:2,:);
 Vnew(3*N/2+1:2*N,:) = -Vbar(1:2*N-3*N/2,:);

 % synthesis
 vv = zeros(1,2*M);
 v = zeros(1,M);
 u = zeros(1,N*L);
 for m=1:L,
   % insert new input block
   vv = [Vnew(:,m).',vv(1:2*M-2*N)];

   % extract valid sub-blocks
   for k=0:M/2/N-1,
     v(k*2*N+1:(k+1)*2*N) = vv(k*4*N+[1:N,3*N+1:4*N]);
   end;

   % window using filter coeffs
   vh = v.*hw;

   % combine
   for k=1:N,
     u((m-1)*N+k) = sum(vh([k:N:M]));
   end;
 end;
