function x=fasthadamard(x)
% This function implement the Hadamard Transform. 
% Code was based on fwht1d.m, written by Gylson Thomas,
% gylson_thomas@yahoo.com, in 2005-2007.
% Code was extensively modified by Matthew Bush, mbbush@gmail.com, in 2013
% to make it parallelizable. 
n = floor(log2(length(x)));
N = pow2(n);
x = x(1:N);
% indexes(1:N/2*n,1:9)=0; % with troubleshooting data
%indexes(1:N/2*n,1:2)=0;  % without troubleshooting data
index1(1:N/2*n)=0;
index2(1:N/2*n)=0;
% k1=N; k2=1; k3=N/2;
% logfile=fopen('hadamard5.csv','w');
% fprintf(logfile,'k1\t k2\t k3\t i1\t i2\t i3\t L1\t i\t j\t x(i)\t x(j)\n');
i1step=pow2(n-1);
% ops(1:N,1:N*n,1:2)=0;
% working in parallel, generate a list of operations to do in serial later
parfor ind=1:N/2*n

    
    i1=ceil(ind/i1step) % i1 runs from 1 to n, increasing by 1 every N/2 iterations
    k1=pow2(n+1-i1); % k1 runs from N to 2, decreasing by factor 2 each time i1 is incremented
  %  k2=pow2(i1-1); % k2 runs from 1 to N/2, increasing by factor 2 each time i1 is incremented
    k3=pow2(n-i1); % k3 runs from N/2 to 1, decreasing by factor 2 each time i1 is incremented
    i2=ceil((mod(ind-1,i1step)+1)/k3); % i2 runs from 1 to k2 increasing by 1 every k3 iterations
                            % and resetting to 1 each time i1 is incremented
    i3=mod(ind-1,k3)+1; % i3 runs from 
    L1=i2*k1-k1;
    i=i3+L1; 
    j=i+k3;
%     indexes(ind,3:9)=[k1,k2,k3,i1,i2,i3,L1]; % troubleshooting data
    index1(ind)=i;index2(ind)=j;
%     ops(i,ind)=j;   % later on, add j to i.
%     ops(j,ind)=-i;
end

% working serially, apply the hadamard transform to each pair of indexes
% specified in the array indexes.
for k=1:length(index1)
           temp1= x(index1(k)); 
%             fprintf(logfile,'%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\t%g\n',...
%                 indexes(k,3),indexes(k,4),indexes(k,5),indexes(k,6),...
%                 indexes(k,7),indexes(k,8),indexes(k,9),...
%                 indexes(k,1),indexes(k,2),x(indexes(k,1)),x(indexes(k,2)));
           x(index1(k)) = temp1 + x(index2(k));
           x(index2(k)) = temp1 - x(index2(k));
end
x=x/sqrt(N); 

