Artyom Grigoryan's Official Site

We first consider the paired fast Fourier transform of order 8. The cases of real-input and complex-input are considered separately.
When an input is real only 2 real multiplications are used, and 4 real multiplications when an input is complex.
% ------------------------------------------------------------------
% call: demo_pfft8.m
% Demostration of the fast paired algorithm for the 8-point DFT with
% a minimum number of real multiplications by the factor of a=0.7071:
% 2 multiplications when an input is real
% 4 multiplications when an input is complex
%
% A.M. Grigoryan, 12/8/2008
%
%1. Real-signal composition and transformation
x=[1,4,2,3,7,2,1,4]';
F=pfft8_inputreal(x);
F1=reshape(F,8,1)
24.0000
-3.8787 - 1.7071i
5.0000 + 1.0000i
-8.1213 + 0.2929i
-2.0000
-8.1213 - 0.2929i
5.0000 - 1.0000i
-3.8787 + 1.7071i
%
%2. Complex-signal composition and transformation
y=[1+2j,4+j,2+8j,3+j,7+4j,2+2j,1+3j,4+j]';
G=pfft8_inputcomplex(y);
G1=reshape(G,8,1);
24.0000 -22.0000i
-8.1716 + 1.0000i
4.0000 + 6.0000i
-2.4142 + 1.5858i
-2.0000 -12.0000i
-13.8284 + 1.0000i
6.0000 + 4.0000i
0.4142 + 4.4142i
% -------------------------------------------------------
% call: pfft8_inputreal.m
% The 8-point forward DFT with 1 operation of multiplication
%
% A.M. Grigoryan, 12/8/2008
%
function F=pfft8_inputreal(x)
% Direct calculation of the paired transform
% Matrix of the 8-point paired transform:
% X8= [1 0 0 0 -1 0 0 0
% 0 1 0 0 0 -1 0 0
% 0 0 1 0 0 0 -1 0
% 0 0 0 1 0 0 0 -1
% 1 0 -1 0 1 0 -1 0
% 0 1 0 -1 0 1 0 -1
% 1 -1 1 -1 1 -1 1 -1
% 1 1 1 1 1 1 1 1];
% y_1=X8*x;
% Fast calculation of the paired transform
y_1=fst_1d(x);
%
F(5)=y_1(7); % [p=4]
F(1)=y_1(8); % [p=0]
y_2=y_1(1:6)'.*[1,1-j,-j,-1-j,1,-j];
F(3)=y_2(5)+y_2(6); % [p=2]
F(7)=conj(F(3)); % [p=6]
a=sqrt(2)/2; % nontrivial factor
g1=y_2(1)+y_2(3);
g2=(y_2(2)+y_2(4))*a; % (!) 2 real multiplications are here
F(6)=g1-g2; % [p=5]
F(2)=g1+g2; % [p=1]
F(4)=conj(F(6)); % [p=3]
F(8)=conj(F(2)); % [p=7]
% -------------------------------------------------------
% call: pfft8_inputcomplex.m
% The 8-point forward DFT with 2 operations of multiplication
%
% A.M. Grigoryan, 12/8/2008
%
function F=pfft8_inputcomplex(x)
a=sqrt(2)/2;
% Direct calculation of the paired transform
% Matrix of the 8-point paired transform:
% X8= [1 0 0 0 -1 0 0 0
% 0 1 0 0 0 -1 0 0
% 0 0 1 0 0 0 -1 0
% 0 0 0 1 0 0 0 -1
% 1 0 -1 0 1 0 -1 0
% 0 1 0 -1 0 1 0 -1
% 1 -1 1 -1 1 -1 1 -1
% 1 1 1 1 1 1 1 1];
% y_1=X8*x;
% Fast calculation of the paired transform
y_1=fst_1d(x);
%
F(5)=y_1(7); % [p=4]
F(1)=y_1(8); % [p=0]
y_2=conj(y_1(1:6)').*[1,1-j,-j,-1-j,1,-j];
F(3)=y_2(5)+y_2(6); % [p=2]
F(7)=y_2(5)-y_2(6); % [p=6]
X4=[1 0 -1 0
0 1 0 -1
1 0 1 0
0 1 0 1];
y_3=X4*conj(y_2(1:4)');
y_3(2)=y_3(2)*a; (!) 2 real multiplications are here
y_3(4)=y_3(4)*a; (!) 2 real multiplications are here
T=[1 j 0 0
1 -j 0 0
0 0 1 -1
0 0 1 1];
nn=[8,4,6,2]; % [p=7,3,5,1]
F(nn)=T*y_3;
% -------------------------------------------------------

Below is the version of the first code for the fast paired FFT.
% call: fft_1d.m
% 1-D fast direct Fourier transform of order N=2**M, M>1.
%
% This program can be used also for the Hadamard transform,
% for what we should elluminate all multiplication by
% twiddle factors W in (*) and omit Step 1.
% code fst_1d.m is used for the fast paired transform
%
% A.M. Grigoryan, March 15-16, 1999
%
function B1=fft_1d(A1,ND)
%1. N/2 Exponentional coefficients for diagonal matrices [W]
N2=ND/2;
W=zeros(1,N2);
T=+j*pi/N2;
J=0;
for I=0:N2-1
W(I+1)=exp(J); J=J+T;
end;
%1.1 The case W=ones(1,N2) results in the Hadamard transform
%2. The decomposition by means of the paired transforms (fst_1d)
MD=log2(ND);
NK1=ND; MD1=MD;
B1=A1;
IW=1;
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); % fast paired transform
LK1=LK; if MD1>1 % this 'if' can be moved
II=IW; JJ=J1;
for I2=1:MD1-1
J2=1;
for J3=1:LK1
B1(JJ)=B1(JJ)*W(J2); % (*)
JJ=JJ+1;
J2=J2+II;
end;
LK1=LK1/2;
II=II*2;
end; end; % can be moved
end;
NK1=LK;
MD1=MD1-1;
IW=IW*2;
end;

Another code for the fast paired transform, which is written in recursive form:
% call: paired_1dfft.m
% M.M. Grigoryan, March 2008
function h=paired_1dfft(x)
N=length(x);
% 1. Calculation of the exponential coefficients. The calculations
% can be performed outside the body of this functions. For instance,
% all coefficients wn can be considered as static variables in the code.
% Note, if all wn=1, then this code results in the Hadamard transform.
m=bitshift(N,-1);
t=0:m-1;
wn=exp(-(pi*j/m)*t);
w=[];
while m>1
w=[w wn];
wn=wn(1:2:end);
m=bitshift(m,-1);
end
w=[w 1 1];
% 2. Calculation of the short DFTs of the modified splitting signals
if N==1
h=x;
else
z=fastpaired_1d(x);
z=z.*w;
nn=1;
nk0=bitshift(N,-1);
nk=nk0;
p=[];
while nk0>1
t=nn:nk;
y=z(t);
nk0=bitshift(nk0,-1);
p=[p paired_1dfft(y)];
nn=nk+1;
nk=nk+nk0;
end
h=[p -z(end-1) z(end)];
end;