function [hn] = poly2(cr)

% Assumes final polynomial has all real coefficients.
% Example: poly2([1 2 3 4 3 2 1 -1 0 -1 0 ]) which returns
% 1   -14    75  -180   123   258  -475    80   276  -144     0     0
% The value at the first index corresponds to the coefficient on the x^11
% term.  So, the constant coefficient and the x coefficient are zero.
 
%%% Copyright (c) 1997-1998 Shawn McCaslin
%%% All Rights Reserved.
%%%
%%% This program is free software; you can redistribute it and/or modify
%%% it under the terms of the GNU General Public License as published by
%%% the Free Software Foundation; either version 2 of the License, or
%%% (at your option) any later version.
%%%
%%% This program is distributed in the hope that it will be useful,
%%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%%% GNU General Public License for more details.
%%%
%%% The GNU Public License is available in the file LICENSE, or you
%%% can write to the Free Software Foundation, Inc., 59 Temple Place -
%%% Suite 330, Boston, MA 02111-1307, USA, or you can find it on the
%%% World Wide Web at http://www.fsf.org.

r = cplxpair(cr);

% Combine complex-conj roots: 

j=1;
k=1;
while(j<=length(cr))
	if(imag(cr(j))~=0)
		cc(k,1:3) = conv([1 -cr(j)],[1 -cr(j+1)]);
		j=j+2;
		k=k+1;
	else
		cc(k,1:3) = [ 0 1 -cr(j)];
		j = j+1;
		k=k+1;
	end;
end;
if(rem(size(cc),2)~=0) cc(k,1:3) = [ 0 0 1 ]; end;

% Assume final polynomial is real:

cc = real(cc);

% Sort by maximum magnitude of frequency response:

for j=1:size(cc)
	mfz(j) = max(freqz(cc(j,1:3)));
end;

[mfzs, I1] = sort(-mfz);

% Combine pairs of roots:

m=0;
while(1)
	m=m+1;

	ct1 = cc(I1(1),1:3);

	I1 = I1(2:length(I1));

	if(length(I1) == 1)
		cq(m,1:5) = conv(ct1,cc(I1(1),1:3));
		break;
	end;

	k=0;
	for j=1:length(I1)
		ct2 = conv(ct1,cc(I1(j),1:3));
		H = freqz(ct2,1,128);
		k(j) = max(abs(H));
	end;

	[minv,mini] = min(k);

	cq(m,1:5) = conv(ct1,cc(I1(mini),1:3));

	if(mini==1)
		I1 = I1(2:length(I1));
	elseif(mini==length(I1))
		I1 = I1(1:length(I1)-1);
	else
		I1 = I1([1:mini-1 mini+1:length(I1)]); 
	end;

	if(length(I1) == 1)
		cq(m+1,1:5) = [0 0 0 0 cc(I1(1),1:5)];
		break;
	end;
end;

% Expand polynomial:

c=1;
while(size(cq)>1)

	k=0;
	for j=1:size(cq)

		H = freqz(conv(c,cq(j,1:5)),1,512);
		K2 = abs(H);
		I1 = find(K2~=0);
		K2 = log10(K2(I1)/max(K2(I1)));
		k(j) = cov(K2-mean(K2));
	end;

	[minv,mini] = min(k);

	c  = conv(c,cq(mini,1:5));

	if(mini==1)
		cq = cq(2:size(cq),1:5);
	elseif(mini==size(cq))
		cq = cq(1:size(cq)-1,1:5);
	else
		cq(mini:size(cq)-1,1:5) = cq(mini+1:size(cq),1:5);
		cq = cq(1:size(cq)-1,1:5);
	end;
end;

c  = conv(c,cq(1,1:5));

% Strip off extraneous zeros:

s = find(c~=0);
hn = c(s(1):length(c));

return;
