function X = DFTPolyAnalysisReal(x,h,K,N);
% [X,x_out] = DFTPolyAnalysisReal(x,h,K,N,x_init);
%
% Oversampled DFT filter bank to produce K/2 complex subbands 
% decimated by a factor N. The real input signal is decomposed 
% by complex bandpass filters derived from a prototype filter h. 
% The length of h must be a multiple of 2M, where M is the least
% common multiple of K and N.
%
% This implementation uses real prototype polyphase components,
% followed by a GDFT modulation.
%
% Input parameters:
%      x      input signal (real)
%             (length must be multiple of N)
%      h      prototype filter (real) 
%      K      number of channels (0 ... 2pi)
%      N      decimation ratio
%      x_in   initialization of state vector
%             (if not given, x_in is zero-padded)
%
% Output parameters:
%      X      rows contain decimated subband signals
%             (complex valued)
%      x_out  last Lh-1 samples of input signal
%             (may be used to initialize next block
%              of data)
%
% initial state conditions not implemented yet!
% St.Weiss, University of Strathclyde, 18.8.1997

% ***** calculate parameters *****
M = lcm(K,N);
L = M/N; J = M/K;
B = N/J;               % # of blocks

% ***** validate input parameters *****
[m,n] = size(h);
if m>n,	               % ensure h is a row vector
  h = h.';
end;
Lh = length(h);
if mod(Lh,2*M) ~= 0,
  error(sprintf('length of prototype filter not a multiple of %d',2*M));
end;  
Lx = length(x);        % check length of input sequence
if mod(Lx,N) ~= 0,
  error('length of input sequence muts be multiple of decimation ratio N');
end;

% ***** polyphase decomposition of input signal *****
[m,n] = size(x);
if n > m,              % ensure x is a column vector
  x = x.';
end;
X = zeros(N,Lx/N);
X = ReArrange(X,x);
X = flipud(X);

% ***** create real valued prototype polyphase components *****
V = zeros(2*M,Lh/2/M);             % create polyphase components in rows
V = ReArrange(V,h.');
V = Expand(V,2*L,'full');        % oversampled by factor 2*L  
for m = 1:2*M,                   % introduce appropriate delay
  V(m,:) = Shift(V(m,:),floor((m-1)/N));
end;  

[Vindex] = DFTPolyEntries2(2*K,N);
Out = zeros(2*K,Lx/N);
for b = 1:B,
  Vindexb = Vindex+b;
  IndexIn = (b:B:N);        % partitioning of polyphase input, J entries
  IndexOut = (b:B:2*K);       % partitioning of polyphase output, L entries 
  for l = 1:2*L,
    for j = 1:J,
      Out(IndexOut(l),:) = Out(IndexOut(l),:) + ...
	  filter(V(Vindexb(l,j),:),1,X(IndexIn(j),:));
    end;  
  end;  
end;   

% ***** GDFT modulation of polyphase filter outputs *****
X = exp(sqrt(-1)*2*(0.5:K/2-.5)'*pi/K*((0.5:2*K-.5)-Lh/2))*Out;