We here consider the Hadamard transformation of order equal to a power of 2. Below are the codes for calculating the fast discrete Hadamard transform (DHdT) by using the paired transform, as well as the traditional algorithm. These codes are illustrated by the following program with an example for the N=8 case.

% --------------------------------------------------------------------
% Call: run_dhdt.m
% Demo code for computing the discrete Hadamard transform (DHdT)
% for orders N equal powers of two.
% List of codes used in this program:
% fht_1d.m - computing the DHdT by the paired transform
% fat_1d.m - computing the DHdT by the traditional method
% had_Matp.m - computing the matrix of the DHdT by paired transform
% had_Mat2.m - computing the matrix of the DHdT
% A.M. Grigoryan, October 1990.
%
N=8;
%1. Integer 8-point test-signal
x=[18,32,3,15,4,31,1,25]';
% Apply both algorithms for the Hadamard transform:
h1=fht_1d(x',N);
h2=fat_1d(x,N);
% the results h1=h2'
% h1=[1,23,25,7,-5,41,-77,129]
%
%2. Calculate the matrix (NxN) of the DHdT
T1=had_Matp(N);
T2=had_Mat2(N);
% the matrices T1=T2 and equal H(N:-1:1,:)
% 1 -1 -1 1 -1 1 1 -1
% 1 1 -1 -1 -1 -1 1 1
% 1 -1 1 -1 -1 1 -1 1
% 1 1 1 1 -1 -1 -1 -1
% 1 -1 -1 1 1 -1 -1 1
% 1 1 -1 -1 1 1 -1 -1
% 1 -1 1 -1 1 -1 1 -1
% 1 1 1 1 1 1 1 1
%
% when comparing with the matrix H of the DHdT from MATLAB
H=hadamard(N);
% 1 1 1 1 1 1 1 1
% 1 -1 1 -1 1 -1 1 -1
% 1 1 -1 -1 1 1 -1 -1
% 1 -1 -1 1 1 -1 -1 1
% 1 1 1 1 -1 -1 -1 -1
% 1 -1 1 -1 -1 1 -1 1
% 1 1 -1 -1 -1 -1 1 1
% 1 -1 -1 1 -1 1 1 -1
y=H*x; % the result y(N:-1:1)=h2=h1'
% y'=[129,-77,41,-5,7,25,23,1]
% ------------------------------------------------ end of the code ----

function B1=fht_1d(A1,ND)
MD=log2(ND);
NK1=ND;
B1=A1;
for I=1:MD
LK=NK1/2;
for J1=1:2*NK1:ND
N_1=NK1+J1-1;
B2=zeros(1,NK1);
B2=B1(J1:N_1);
B1(J1:N_1)=fst_1d(B2);
end;
NK1=LK;
end;
% ------------------------------------------------
function T=had_Matp(N)
T=zeros(N);
for i1=1:N
y=zeros(1,N);
y(i1)=1;
T(:,i1)=fht_1d(y,N);
end;
% ------------------------------------------------
function B1=fat_1d(A1,ND)
MD=log2(ND);
IK=1; if(IK == 1) D=1.; else D=2.; end;
B1=A1;
NK1=ND;
for I=1:MD
LK=NK1/2;
for J=1:NK1:ND
N_1=LK+J-1;
for K=J:N_1
K1=K+LK; T2=B1(K); T1=B1(K1);
B1(K)=(T2-T1)/D; B1(K1)=(T2+T1)/D;
end;
end;
NK1=LK; IK=-IK;
end;
% ------------------------------------------------
function T=had_Mat2(N)
T=zeros(N);
for i1=1:N
y=zeros(1,N);
y(i1)=1;
T(:,i1)=fat_1d(y,N);
end;
% ------------------------------------------------

(!) This is a very interesting fact which shows a strong relation between the Hadamard and Fourier transforms. The fast calculation of the Hadamard trasform (in the N=2^r case) can be obtained from the fast Fourier transform. Here we consider the fast algorithms of the Hadamard transform obtained from the paired fast Fourier transform. The method is very simple: all rotation coefficients W^k in the paired algorithm are considered to be 1. Below is the demo program wherein two codes of the paired transform-based FFT are used after simplification of these codes for the case when N is a power of two.

% Call: run_dhdt_byfourier.m
% Demo code for computing the matrix of the Hadamard transform
% by simplifying codes of the paired transform-based fast Fouirer transform.
% List of codes used here:
% fht_1d.m - computing the DHdT by the paired-Foureir-transform
% paired_1dFH.m - computing the DHdT by the paired-Foureir-transform
% dhd_matrix1.m - computing the matrix of the DHdT by 'fht_1d.m'
% dhd_matrix2.m - computing the matrix of the DHdT by 'paired_1dFH.m'
% Note that "fht_1d.m" uses the original fast paired transform, as well as "paired_1dFH.m" which has a recurrent structure.
% Artyom M. Grigoryan, December 2008.

N=8;
HM8a=dhd_matrix1(N);
1 -1 -1 1 -1 1 1 -1
1 1 -1 -1 -1 -1 1 1
1 -1 1 -1 -1 1 -1 1
1 1 1 1 -1 -1 -1 -1
1 -1 -1 1 1 -1 -1 1
1 1 -1 -1 1 1 -1 -1
1 -1 1 -1 1 -1 1 -1
1 1 1 1 1 1 1 1

HM8b=dhd_matrix2(N);
% In this case we have the complex Hadamard matrix:
[1 i -1 -i -1 -i 1 i
1 -i -1 i -1 i 1 -i
1 -1 1 -1 -1 1 -1 1
1 1 1 1 -1 -1 -1 -1
1 -1 -1 1 1 -1 -1 1
1 1 -1 -1 1 1 -1 -1
1 -1 1 -1 1 -1 1 -1
1 1 1 1 1 1 1 1];

% We can verify if this complex matrix is a unitary Hadamard matrix:
HM8b*HM8b'
8 0 0 0 0 0 0 0
0 8 0 0 0 0 0 0
0 0 8 0 0 0 0 0
0 0 0 8 0 0 0 0
0 0 0 0 8 0 0 0
0 0 0 0 0 8 0 0
0 0 0 0 0 0 8 0
0 0 0 0 0 0 0 8
% Yes it is.
% ----------------------------------------------------------------