% Polyphase Quadrature Filterbank: Analysis
%
% usage:
%   S = pqf_anal(x,h,N);
% where:
%   N	  : number of subband branches
%   h     : coefficients of lowpass prototype filter (cutoff at pi/2N radians)
%   x     : vector of time-domain input samples
%   S     : matrix of decimated subband outputs (having N rows)
%   
% Copyright 1999 --- Phil Schniter   

function S = pqf_anal(x,h,N);

 % determine constants
 M = length(h);
 L = length(x);
 if N*floor(M/N)~=M, error('Filter length must be a multiple of N!'); end;

 % zero-pad input for convenience
 xp = [zeros(1,N-1),x(:).',zeros(1,N-mod(L+N-1,N)),zeros(1,M)];

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

 % analysis
 W = zeros(2*N,ceil(length(xp)/N));
 xx = zeros(1,M);
 for m=1:ceil(length(xp)/N),
   % load new input block
   xx = [xp((m-1)*N+[N:-1:1]), xx(1:M-N)];

   % window using filter coeffs
   xh = xx.*hw;

   % combine
   for k=1:2*N,
     W(k,m) = sum(xh(k:2*N:M));
   end;
 end; 

 % IDCT
 Wbar = zeros(N,ceil(length(xp)/N));
 Wbar(1,:) = sqrt(2)*W(N/2+1,:);
 Wbar(2:N/2+1,:) = W(N/2+2:2*N/2+1,:) + W(N/2:-1:1,:);
 Wbar(N/2+2:N,:) = W(N+2:N/2+N,:) - W(2*N:-1:N+N/2+2,:);
 S = sqrt(N/2)*idct(Wbar);
