function [B,A] = warppoles(a,alpha)
%  [B,A] = warppoles(a,alpha)  warp an all-pole polynomial by substitution
%    Warp all-pole polynomials defined by rows of a by the first-order 
%    warp factor alpha.  Negative alpha shifts poles up in frequency.
%    Output polynomials have zeros too, hence B and A.
% 2003-12-10 dpwe@ee.columbia.edu

% Construct z-hat^-1 polynomial
d = [-alpha 1];
c = [1 -alpha];

[nrows,order] = size(a);

A = zeros(nrows, order);
B = zeros(nrows, order);

B(:,1) = a(:,1);
A(:,1) = ones(nrows,1);

dd = d;
cc = c;

% This code originally mapped zeros.  I adapted it to map
% poles just by interchanging b and a, then swapping again at the 
% end.  Sorry that makes the variables confusing to read.

for n = 2:order

  for row = 1:nrows

    % add another factor to num, den
    B(row,1:order) = conv(B(row,1:(order-1)), c);

  end

  % accumulate terms from this factor
  B(:,1:length(dd)) = B(:,1:length(dd)) + a(:,n)*dd;
    
  dd = conv(dd, d);
  cc = conv(cc, c);

end
    
% Construct the uniform A polynomial (same for all rows)
AA = 1;
for n = 2:order
  AA = conv(AA,c);
end
A = repmat(AA, nrows, 1);
  

% Exchange zeros and poles
T = A; A = B; B = T;