Artyom Grigoryan's Official Site

The universal fast algorithm, for computing the 1-D discrete Fourier transform (DFT) and discrete Hadamard transform (DHdT), is presented below. The order of the transforms is the power of two. The algorithm is based on the concept of the paired transform which transfers the 1-D signal into the frequency-and-time domain. The signal is represented by a set of log2(N) short signals generated by frequencies 1,2,4,8,..., N/2, and 0. The length of these splitting-signals are N/2,N/4,N/8,N/16,...,1, and 1, respectively. By my opinion, this is the best algorithm for implementating calculations of the 1-D DFT and DHdT on the computer or FPGA. / Art Grigoryan

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 multiplications 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;