function y = DFTPolySynthesisCmplx(Y,h,K,N);
% [y,Y_out] = DFTPolySynthesisCmplx(Y,h,K,N,Y_init);
%
% Oversampled DFT synthesis filter bank to reconstruct a real 
% signal from K complex valued subbands decimated by a factor N. 
% The complex bandpass filters are 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.
%
% The subband signals are first rotated by a GDFT matrix. Then, the
% real values are processed by real prototype polyphase filters. 
%
% Input parameters:
%      Y      K complex subband signals in rows
%      h      prototype filter (real) 
%      K      number of channels (0 ... 2pi)
%      N      decimation ratio
%      Y_init initialization of state vectors
%             (if not given, Y_init is zero-padded)
%
% Output parameters:
%      y      reconstructed real fullband signal
%      Y_out  last Lh/N-1 samples of subband signals
%             (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 = size(Y,2)*N;

% ***** take real part of GDFT modulation of subband signals *****
Z = exp(sqrt(-1)*2*(0.5:K-.5)'*pi/K*((0.5:2*K-.5)-Lh/2))'*Y;

% ***** create polyphase components of real valued prototype filter *****
V = zeros(2*M,Lh/2/M);           % create polyphase components in rows
V = ReArrange(V,h.'/sqrt(2));
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);
Vt = fliplr(conj(V));

Y2 = zeros(N,Lx/N);
for b = 1:B,
  Vindexb = Vindex+b;
  IndexIn = (b:B:2*K);           % L entries
  IndexOut = (b:B:N);            % J entries
  for j = 1:J,
    for l = 1:2*L,
      Y2(IndexOut(j),:) = Y2(IndexOut(j),:) + ...
	  filter(Vt(Vindexb(l,j),:),1,Z(IndexIn(l),:));
    end;
  end;
end;
Y2 = flipud(Y2);
y = 2*N*Y2(:);