function A=halfgen4(up,N);
% up is the stopband width, as a fraction of input sampling rate
% N is the order of half-band filter to generate
% A is the full set of FIR coefficients, 4*N-1 long

npt=N*20;

wmax=2*pi*up;
x0=([0:npt]-.0)'/npt;
yfit=1-x0.^2;   % possibly bogus, but good enough to get started
wfit=yfit*wmax;

q=[1:2:(2*N-1)];
target=.5*ones(length(wfit),1);
basis=cos(wfit*q);
l=basis\target;

% gset style data linespoints
% hold off
% plot(wfit,basis*l-.5);

% Not Pretty.  I got frustrated and just brute-forced the weights.
% I have no reason to think the process is stable.
weight=1+wfit*0;
hold on
for iter=[1:35]
    err=abs(basis*l-.5);
    maxerr=max(err);
    ix=find(err>maxerr*0.99);
    boost=1+1.5/(iter+10);
    weight(ix)=weight(ix)*boost;

    wbasis=basis.*(weight*ones(1,N));
    wtarget=target.*weight;
    l=wbasis\wtarget;
    % plot(wfit,basis*l-.5);
end

err=abs(basis*l-.5);
maxerr=max(err);
ix=find(err>maxerr*0.99);
%hold off
%plot(wfit,abs(basis*l-.5),"-",wfit(ix),abs(basis(ix,:)*l-.5),"*",[0;wmax],[maxerr,maxerr],"-");

% Expand the N-long l array to the (4*N-1)-long A array,
% suitable for polyval(A,z) computation of filter behavior
A=reshape([l';l'*0],1,length(l)*2);
A=A(1:length(A)-1);
A=[fliplr(A) 1 A]/2;