function xt=fwht(x,blksiz)

% xt = FWHT (x,blksiz) recursively calculates the Walsh-Hadamard transform
% with results in natural order
%
% x:  data to transform, as a row vector
% blksiz: block size for x (defaults to 1)
%
% length(x)/blksiz is 2^(order of the transform)
% elements x(1), x(2),..., x(blksiz) taken as first block
% elements x(blksiz+1), ..., x(2*blksiz) taken as second block
%
% output explained (xb(j)=jth block of xb):
% first block of xt:   xb(1)+xb(2)+xb(3)+xb(4)...
% second block of xt:  first half of x - second half of x
% third block of xt:   first 1/4 of x - second 1/4 of x + third 1/4 of x - fourth 1/4 of x
%
error(nargchk(1,2,nargin));
if (nargin == 1)
   blksiz=1;
end
lenx=length(x);
blocks=lenx/blksiz;
curord=floor(log2(blocks));
if (curord <= 0)
   xt=x;
   return
else
   xr=transpose(reshape(x,blksiz,blocks)); %routine to cut into even and odd blocks
   x1r=xr(linspace(1,blocks-1,blocks/2),:); 
   x2r=xr(linspace(2,blocks,blocks/2),:);
   x1=reshape(transpose(x1r),1,lenx/2);
   x2=reshape(transpose(x2r),1,lenx/2);
   xt1=fwht(x1,blksiz);
   xt2=fwht(x2,blksiz);
   xt=[xt1+xt2 xt1-xt2];
end