function [Y,H] = dcells(X,zpf,W,aliasbands,debug)
% DCELLS - Break up complex spectrum X into dyadic cells

% Origin: http://ccrma.stanford.edu/~jos/sasp/dcells.m
%
% Reference:
% http://ccrma.stanford.edu/~jos/sasp/Audio_Filter_Banks.html
%
% Example Usage:
%   http://ccrma.stanford.edu/~jos/sasp/Spectral_Rotation_Real_Signals.html

if nargin<5, debug=0; end
if nargin<4, aliasbands=1; end
if nargin<3, W=1; end % (1 means "impulse in freq domain")
if nargin<2, zpf=1; end
if zpf>1 && W==1
  error('dcells: Must compute asinc kernel for given zero-padding factor');
end
X = X(:);
W = W(:)/sum(W); 
N0zp = length(X);
N0 = N0zp/zpf;
NS = round(log2(N0)); % Number of octaves
Y = cell(1,NS+1);     % extra cell for lowest band
H = zeros(NS+1,N0zp); % collection of filters (spectral windows)
N = N0 * 2 .^ -(1:NS+1); % Spectral widths in (critically sampled) bins
N(NS+1)=N(NS);        % dc channel looks like next band up
Nzp = N*zpf;
for k=1:NS+1
  if k==NS+1, Hideal = []; else Hideal = zeros(Nzp(k),1); end
  Hideal = [Hideal;ones(Nzp(k),1)];
  if k>1, Hideal = [Hideal;zeros(N0zp-length(Hideal),1)]; end
  if W==1 && zpf==1 && aliasbands!=0
    % no overlap-add decomposition kernel specified, so
    % simply partition the FFT bins dyadically:
    if k == (NS+1)
      Y{k} = X(1); % dc band
    else
      Y{k} = X(N(k)+1:2*N(k)) .'; % top half spectrum = top octave
    end
    H(k,:) = Hideal;
  else
    % Need to build an overlap-add spectral window.
    % Note that these filters H can be precomputed:
    H(k,:) = cconvr(W,Hideal,N0zp); % real circular convolution
    % cconvr() available as http://ccrma.stanford.edu/~jos/sasp/cconvr.m
    toshow = [W/W(1),Hideal,H(k,:).'];
    if debug>1
      plot(toshow,'-*'); legend('Kernel','Hideal','H'); 
      title(sprintf('Filter for Band %d (linear)',k)); grid('on');
      pause;
      plot(db(toshow,-100),'-*'); grid('on'); legend('Kernel','Hideal','H'); ...
      title(sprintf('Filter for Band %d (dB)',k)); grid('on');
      pause;
    end
    Yek = (X.') .* H(k,:); % Apply band filter as a spectral window
                           % disp(sprintf("N(k) = %d",N(k)));
    if aliasbands
      Y{k} = alias(Yek,N0zp/Nzp(k)); % Aliasing onto a block of length Nzp(k)
      % alias() available as http://ccrma.stanford.edu/~jos/sasp/alias.m
    else
      Y{k} = Yek; % Filter bank retains perfect reconstruction property
    end
    if debug>2
      plot(abs(Y{k}),'-*'); 
      title(sprintf('|Spectral Band| %d',k)); grid('on');
      pause; 
    end    
  end
end
Hsum = sum(H);
fbserr = norm(Hsum-ones(size(Hsum)));
disp(sprintf('Filter bank summation L2 error = %0.2f',fbserr));
if debug && fbserr > 1E-7
  plot(Hsum,'-*'); title('Filter Bank Sum'); 
  disp('PAUSING'); pause;
end
if debug>1
  plot(H','-*',"linewidth",10); 
  title('Complex Octave Filter Bank',"fontsize",18);
  axis([1, length(H), 0, max(max(abs(H)))*1.1]);
  bandslegend(NS+1,"'east'");
  xlabel('Spectral Sample Number',"fontsize",18);
  ylabel('Filter Gain',"fontsize",18);
  disp('PAUSING'); pause;
end
